]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
new FE_DGP
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 11 Feb 2002 12:58:57 +0000 (12:58 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 11 Feb 2002 12:58:57 +0000 (12:58 +0000)
git-svn-id: https://svn.dealii.org/trunk@5501 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/fe/fe_dgp.h [new file with mode: 0644]
deal.II/deal.II/source/fe/fe_dgp.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/fe/fe_dgp.h b/deal.II/deal.II/include/fe/fe_dgp.h
new file mode 100644 (file)
index 0000000..787ea5f
--- /dev/null
@@ -0,0 +1,364 @@
+//---------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 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_dgp_h
+#define __deal2__fe_dgp_h
+
+#include <base/config.h>
+#include <base/polynomial.h>
+#include <fe/fe.h>
+
+template <int dim> class PolynomialSpace;
+template <int dim> class MappingQ;
+
+
+/**
+ * Discontinuous tensor product elements based on equidistant support points.
+ *
+ * This is a discontinuous finite element using interpolating tensor
+ * product polynomials. The shape functions are Lagrangian
+ * interpolants of an equidistant grid of points on the unit cell. The
+ * points are numbered in lexicographical order, @p{x} running fastest.
+ *
+ * @author Guido Kanschat, Ralf Hartmann, 2001
+ */
+template <int dim>
+class FE_DGP : public FiniteElement<dim>
+{
+  public:
+                                    /**
+                                     * Constructor for tensor product
+                                     * polynomials of degree @p{k}.
+                                     */
+    FE_DGP (unsigned int k);
+                                    /**
+                                     * Destructor.
+                                     */
+    ~FE_DGP ();
+    
+                                    /**
+                                     * Return the value of the
+                                     * @p{i}th shape function at the
+                                     * point @p{p}.  @p{p} is a point
+                                     * on the reference element.
+                                     */
+    virtual double shape_value (const unsigned int i,
+                               const Point<dim> &p) const;
+    
+                                    /**
+                                     * Return the gradient of the
+                                     * @p{i}th shape function at the
+                                     * point @p{p}. @p{p} is a point
+                                     * on the reference element, and
+                                     * likewise the gradient is the
+                                     * gradient on the unit cell with
+                                     * respect to unit cell
+                                     * coordinates.
+                                     */
+    virtual Tensor<1,dim> shape_grad (const unsigned int  i,
+                                     const Point<dim>   &p) const;
+
+                                    /**
+                                     * Return the tensor of second
+                                     * derivatives of the @p{i}th
+                                     * shape function at point @p{p}
+                                     * on the unit cell. The
+                                     * derivatives are derivatives on
+                                     * the unit cell with respect to
+                                     * unit cell coordinates.
+                                     */
+    virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
+                                          const Point<dim> &p) const;
+    
+                                    /**
+                                     * Return the polynomial degree
+                                     * of this finite element,
+                                     * i.e. the value passed to the
+                                     * constructor.
+                                     */
+    unsigned int get_degree () const;
+
+                                    /**
+                                     * Number of base elements in a
+                                     * mixed discretization. Since
+                                     * this is a scalar element,
+                                     * return one.
+                                     */
+    virtual unsigned int n_base_elements () const;
+    
+                                    /**
+                                     * Access to base element
+                                     * objects. Since this element is
+                                     * scalar, @p{base_element(0)} is
+                                     * @p{this}, and all other
+                                     * indices throw an error.
+                                     */
+    virtual const FiniteElement<dim> & base_element (const unsigned int index) const;
+    
+                                    /**
+                                     * Check for non-zero values on a face.
+                                     *
+                                     * This function returns
+                                     * @p{true}, if the shape
+                                     * function @p{shape_index} has
+                                     * non-zero values on the face
+                                     * @p{face_index}.
+                                     *
+                                     * Implementation of the
+                                     * interface in
+                                     * @ref{FiniteElement}
+                                     */
+    virtual bool has_support_on_face (const unsigned int shape_index,
+                                     const unsigned int face_index) const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object.
+                                     *
+                                     * This function is made virtual,
+                                     * since finite element objects
+                                     * are usually accessed through
+                                     * pointers to their base class,
+                                     * rather than the class itself.
+                                     */
+    virtual unsigned int memory_consumption () const;
+
+  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
+                                     * @ref{FiniteElement}.
+                                     */
+    virtual void
+    fill_fe_values (const Mapping<dim> &mapping,
+                   const typename DoFHandler<dim>::cell_iterator &cell,
+                   const Quadrature<dim>                &quadrature,
+                   typename Mapping<dim>::InternalDataBase      &mapping_internal,
+                   typename Mapping<dim>::InternalDataBase      &fe_internal,
+                   FEValuesData<dim>& data) const;
+    
+                                    /**
+                                     * Implementation of the same
+                                     * function in
+                                     * @ref{FiniteElement}.
+                                     */
+    virtual void
+    fill_fe_face_values (const Mapping<dim> &mapping,
+                        const typename DoFHandler<dim>::cell_iterator &cell,
+                        const unsigned int                    face_no,
+                        const Quadrature<dim-1>                &quadrature,
+                        typename Mapping<dim>::InternalDataBase      &mapping_internal,
+                        typename Mapping<dim>::InternalDataBase      &fe_internal,
+                        FEValuesData<dim>& data) const ;
+    
+                                    /**
+                                     * Implementation of the same
+                                     * function in
+                                     * @ref{FiniteElement}.
+                                     */
+    virtual void
+    fill_fe_subface_values (const Mapping<dim> &mapping,
+                           const typename DoFHandler<dim>::cell_iterator &cell,
+                           const unsigned int                    face_no,
+                           const unsigned int                    sub_no,
+                           const Quadrature<dim-1>                &quadrature,
+                           typename Mapping<dim>::InternalDataBase      &mapping_internal,
+                           typename Mapping<dim>::InternalDataBase      &fe_internal,
+                           FEValuesData<dim>& data) const ;
+
+  private:
+
+                                    /**
+                                     * Declare a nested class which
+                                     * will has static definitions of
+                                     * various matrices such as
+                                     * constraint and embedding
+                                     * matrices. The definition of
+                                     * the various static fields are
+                                     * in the files @p{fe_q_[123]d.cc}
+                                     * in the source directory.
+                                     */
+    struct Matrices
+    {
+                                        /**
+                                         * Pointers to the embedding
+                                         * matrices, one for each
+                                         * polynomial degree starting
+                                         * from constant elements
+                                         */
+       static const double * const embedding[];
+
+                                        /**
+                                         * Number of elements (first
+                                         * index) the above field
+                                         * has. Equals the highest
+                                         * polynomial degree plus one
+                                         * for which the embedding
+                                         * matrices have been
+                                         * computed.
+                                         */
+       static const unsigned int n_embedding_matrices;
+
+                                        /**
+                                         * As @p{embedding} but for
+                                         * projection matrices.
+                                         */
+       static const double * const projection_matrices[];
+
+                                        /**
+                                         * As
+                                         * @p{n_embedding_matrices}
+                                         * but for projection
+                                         * matrices.
+                                         */
+       static const unsigned int n_projection_matrices;
+    };
+
+    
+                                    /**
+                                     * 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(unsigned int degree);
+    
+                                    /**
+                                     * Compute flags for initial update only.
+                                     */
+    virtual UpdateFlags update_once (const UpdateFlags flags) const;
+  
+                                    /**
+                                     * Compute flags for update on each cell.
+                                     */
+    virtual UpdateFlags update_each (const UpdateFlags flags) const;
+  
+                                    /**
+                                     * Degree of the polynomials.
+                                     */  
+    const unsigned int degree;
+
+                                    /**
+                                     * Pointer to the tensor
+                                     * product polynomials.
+                                     */
+    PolynomialSpace<dim>* poly;
+
+                                    /**
+                                     * Fields of cell-independent data.
+                                     */
+    class InternalData : public FiniteElementBase<dim>::InternalDataBase
+    {
+      public:
+                                        /**
+                                         * Array with shape function
+                                         * values in quadrature
+                                         * points. There is one
+                                         * vector for each shape
+                                         * function, containing
+                                         * values for each quadrature
+                                         * point.
+                                         */
+       std::vector<std::vector<double> > shape_values;
+       
+                                        /**
+                                         * Array with shape function
+                                         * gradients in quadrature
+                                         * points. There is one
+                                         * vector for each shape
+                                         * function, containing
+                                         * values for each quadrature
+                                         * point.
+                                         */                                  
+       typename std::vector<std::vector<Tensor<1,dim> > > shape_gradients;
+    };
+    
+                                    /**
+                                     * Allow access from other dimensions.
+                                     */
+    template <int dim1> friend class FE_DGP;
+
+                                    /**
+                                     * Allows @p{MappingQ} class to
+                                     * access to build_renumbering
+                                     * function.
+                                     */
+    friend class MappingQ<dim>;
+};
+
+
+// declaration of explicit specializations
+
+template <> 
+const double * const FE_DGP<1>::Matrices::embedding[];
+
+template <>
+const unsigned int FE_DGP<1>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_DGP<1>::Matrices::projection_matrices[];
+
+template <>
+const unsigned int FE_DGP<1>::Matrices::n_projection_matrices;
+
+template <> 
+const double * const FE_DGP<2>::Matrices::embedding[];
+
+template <>
+const unsigned int FE_DGP<2>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_DGP<2>::Matrices::projection_matrices[];
+
+template <>
+const unsigned int FE_DGP<2>::Matrices::n_projection_matrices;
+
+template <> 
+const double * const FE_DGP<3>::Matrices::embedding[];
+
+template <>
+const unsigned int FE_DGP<3>::Matrices::n_embedding_matrices;
+
+template <>
+const double * const FE_DGP<3>::Matrices::projection_matrices[];
+
+template <>
+const unsigned int FE_DGP<3>::Matrices::n_projection_matrices;
+
+
+#endif
diff --git a/deal.II/deal.II/source/fe/fe_dgp.cc b/deal.II/deal.II/source/fe/fe_dgp.cc
new file mode 100644 (file)
index 0000000..80372c2
--- /dev/null
@@ -0,0 +1,411 @@
+//----------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002 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.h>
+#include <base/polynomial.h>
+#include <base/polynomial_space.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/mapping_q1.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_values.h>
+
+
+
+template <int dim>
+FE_DGP<dim>::FE_DGP (unsigned int degree)
+               :
+               FiniteElement<dim> (FiniteElementData<dim>(get_dpo_vector(degree),1),
+                                   std::vector<bool>(1,true)),
+               degree(degree),
+               poly(0)
+{
+                                  // create array of Legendre polynomials
+  std::vector<Legendre<double> > v;
+  for (unsigned int i=0;i<=degree;++i)
+    v.push_back(Legendre<double>(i));
+  poly = new PolynomialSpace<dim> (v);
+
+                                  // if defined, copy over matrices
+                                  // from precomputed arrays
+//    if ((degree < Matrices::n_embedding_matrices) &&
+//        (Matrices::embedding[degree] != 0))
+//      {
+//        prolongation[0].fill (Matrices::embedding[degree]);
+//      }
+//    else
+//                                  // matrix undefined, set size to zero
+//      for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
+//        prolongation[i].reinit(0);
+
+//                                // same as above: copy over matrix
+//                                // from predefined values and
+//                                // generate all others by rotation
+//    if ((degree < Matrices::n_projection_matrices) &&
+//        (Matrices::projection_matrices[degree] != 0))
+//      {
+//        restriction[0].fill (Matrices::projection_matrices[degree]);
+//      }
+//    else
+//                                  // matrix undefined, set size to zero
+//      for (unsigned int i=0;i<GeometryInfo<dim>::children_per_cell;++i)
+//        restriction[i].reinit(0);  
+};
+
+
+
+template <int dim>
+FE_DGP<dim>::~FE_DGP ()
+{
+                                  // delete poly member and set it to
+                                  // zero to prevent accidental use
+  delete poly;
+  poly = 0;
+}
+
+
+
+template <int dim>
+FiniteElement<dim> *
+FE_DGP<dim>::clone() const
+{
+  return new FE_DGP<dim>(degree);
+}
+
+
+
+template <int dim>
+double
+FE_DGP<dim>::shape_value (const unsigned int i,
+                         const Point<dim> &p) const
+{
+  return poly->compute_value(i, p);
+}
+
+
+
+template <int dim>
+Tensor<1,dim>
+FE_DGP<dim>::shape_grad (const unsigned int i,
+                        const Point<dim> &p) const
+{
+  return poly->compute_grad(i, p);
+}
+
+
+
+template <int dim>
+Tensor<2,dim>
+FE_DGP<dim>::shape_grad_grad (const unsigned int i,
+                             const Point<dim> &p) const
+{
+  return poly->compute_grad_grad(i, p);
+}
+
+
+//----------------------------------------------------------------------
+// Auxiliary functions
+//----------------------------------------------------------------------
+
+
+template <int dim>
+std::vector<unsigned int>
+FE_DGP<dim>::get_dpo_vector(unsigned int deg)
+{
+  std::vector<unsigned int> dpo(dim+1, static_cast<unsigned int>(0));
+  dpo[dim] = ++deg;
+  for (unsigned int i=1;i<dim;++i)
+    {
+      dpo[dim] *= deg+i;
+      dpo[dim] /= i+1;
+    }
+  return dpo;
+}
+
+
+template <int dim>
+UpdateFlags
+FE_DGP<dim>::update_once (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+
+  if (flags & update_values)
+    out |= update_values;
+
+  return out;
+}
+
+
+template <int dim>
+UpdateFlags
+FE_DGP<dim>::update_each (const UpdateFlags flags) const
+{
+  UpdateFlags out = update_default;
+
+  if (flags & update_gradients)
+    out |= update_gradients | update_covariant_transformation;
+
+  if (flags & update_second_derivatives)
+    out |= update_second_derivatives | update_covariant_transformation;
+
+  return out;
+}
+
+
+//----------------------------------------------------------------------
+// Data field initialization
+//----------------------------------------------------------------------
+
+template <int dim>
+typename Mapping<dim>::InternalDataBase *
+FE_DGP<dim>::get_data (const UpdateFlags      update_flags,
+                      const Mapping<dim>    &mapping,
+                      const Quadrature<dim> &quadrature) const
+{
+  InternalData* data = new InternalData;
+  std::vector<double> values(0);
+  std::vector<Tensor<1,dim> > grads(0);
+  std::vector<Tensor<2,dim> > grad_grads(0);
+  
+                                  // check what needs to be
+                                  // initialized only once and what
+                                  // on every cell/face/subface we
+                                  // visit
+  data->update_once = update_once(update_flags);
+  data->update_each = update_each(update_flags);
+  data->update_flags = data->update_once | data->update_each;
+
+  const UpdateFlags flags(data->update_flags);
+  const unsigned int n_q_points = quadrature.n_quadrature_points;
+  
+  if (flags & update_values)
+    {
+      values.resize (dofs_per_cell);
+      data->shape_values.resize(dofs_per_cell,
+                               std::vector<double>(n_q_points));
+    }
+
+  if (flags & update_gradients)
+    {
+      grads.resize (dofs_per_cell);
+      data->shape_gradients.resize(dofs_per_cell,
+                                  std::vector<Tensor<1,dim> >(n_q_points));
+    }
+
+                                  // if second derivatives through
+                                  // finite differencing is required,
+                                  // then initialize some objects for
+                                  // that
+  if (flags & update_second_derivatives)
+    data->initialize_2nd (this, mapping, quadrature);
+  
+  
+  if (flags & (update_values | update_gradients))
+    for (unsigned int i=0; i<n_q_points; ++i)
+      {
+       poly->compute(quadrature.point(i), values, grads, grad_grads);
+       for (unsigned int k=0; k<dofs_per_cell; ++k)
+         {
+           if (flags & update_values)
+             data->shape_values[k][i] = values[k];
+           if (flags & update_gradients)
+             data->shape_gradients[k][i] = grads[k];
+         }
+      }
+  return data;
+}
+
+
+
+//----------------------------------------------------------------------
+// Fill data of FEValues
+//----------------------------------------------------------------------
+
+template <int dim>
+void
+FE_DGP<dim>::fill_fe_values (const Mapping<dim>                   &mapping,
+                            const typename DoFHandler<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
+                                  // an exception if that is not
+                                  // possible
+  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+  
+  const UpdateFlags flags(fe_data.current_update_flags());
+
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
+    {
+      for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i];
+      
+      if (flags & update_gradients)
+       mapping.transform_covariant(data.shape_gradients[k],
+                                   fe_data.shape_gradients[k],
+                                   mapping_data, 0);
+    }
+  
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, 0, mapping_data, fe_data, data);
+  
+  fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+void
+FE_DGP<dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
+                                 const typename DoFHandler<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
+                                  // an exception if that is not
+                                  // possible
+  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                  // faces are stored contiguously)
+  const unsigned int offset = face * quadrature.n_quadrature_points;
+  
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
+    {
+      for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+      
+      if (flags & update_gradients)
+       mapping.transform_covariant(data.shape_gradients[k],
+                                   fe_data.shape_gradients[k],
+                                   mapping_data, offset);
+    }
+
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+  
+  fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+void
+FE_DGP<dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
+                                    const typename DoFHandler<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
+                                  // an exception if that is not
+                                  // possible
+  InternalData &fe_data = dynamic_cast<InternalData &> (fedata);
+  
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                  // sub-faces are stored contiguously)
+  const unsigned int offset = (face * GeometryInfo<dim>::subfaces_per_face + subface)
+                             * quadrature.n_quadrature_points;
+
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  for (unsigned int k=0; k<dofs_per_cell; ++k)
+    {
+      for (unsigned int i=0;i<quadrature.n_quadrature_points;++i)
+       if (flags & update_values)
+         data.shape_values(k,i) = fe_data.shape_values[k][i+offset];
+      
+      if (flags & update_gradients)
+       mapping.transform_covariant(data.shape_gradients[k],
+                                   fe_data.shape_gradients[k],
+                                   mapping_data, offset);
+    }
+  
+  if (flags & update_second_derivatives)
+    compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+  
+  fe_data.first_cell = false;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_DGP<dim>::n_base_elements () const
+{
+  return 1;
+};
+
+
+
+template <int dim>
+const FiniteElement<dim> &
+FE_DGP<dim>::base_element (const unsigned int index) const
+{
+  Assert (index==0, ExcIndexRange(index, 0, 1));
+  return *this;
+};
+
+
+
+template <int dim>
+bool
+FE_DGP<dim>::has_support_on_face (const unsigned int shape_index,
+                                 const unsigned int face_index) const
+{
+  return true;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_DGP<dim>::memory_consumption () const
+{
+  Assert (false, ExcNotImplemented ());
+  return 0;
+}
+
+
+
+template <int dim>
+unsigned int
+FE_DGP<dim>::get_degree () const
+{
+  return degree;
+};
+
+
+
+
+template class FE_DGP<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.