From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Mon, 11 Feb 2002 12:58:57 +0000 (+0000)
Subject: new FE_DGP
X-Git-Tag: v8.0.0~18340
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3fa9183a80ea97961f23993eef26fda801a2b8b5;p=dealii.git

new FE_DGP


git-svn-id: https://svn.dealii.org/trunk@5501 0785d39b-7218-0410-832d-ea1e28bc413d
---

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
index 0000000000..787ea5fb30
--- /dev/null
+++ b/deal.II/deal.II/include/fe/fe_dgp.h
@@ -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
index 0000000000..80372c2c5c
--- /dev/null
+++ b/deal.II/deal.II/source/fe/fe_dgp.cc
@@ -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>;