--- /dev/null
+//---------------------------------------------------------------
+// $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
--- /dev/null
+//----------------------------------------------------------------
+// $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>;
projects / dealii.git / commitdiff