]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
checking in preliminary version
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 24 May 2005 06:51:17 +0000 (06:51 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 24 May 2005 06:51:17 +0000 (06:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@10717 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/deal.II/include/fe/fe_poly_tensor.h b/deal.II/deal.II/include/fe/fe_poly_tensor.h
new file mode 100644 (file)
index 0000000..f18a006
--- /dev/null
@@ -0,0 +1,312 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 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_poly_tensor_h
+#define __deal2__fe_poly_tensor_h
+
+
+#include <lac/full_matrix.h>
+#include <fe/fe.h>
+
+/*!@addtogroup fe */
+/*@{*/
+
+/**
+ * This class gives a unified framework for the implementation of
+ * FiniteElement classes based on Tensor valued polynomial spaces like
+ * PolynomialsBDM and PolynomialsRaviartThomas.
+ *
+ * Every class that implements following function can be used as
+ * template parameter POLY.
+ *
+ * @code
+ * void compute (const Point<dim>            &unit_point,
+ *               std::vector<Tensor<1,dim> > &values,
+ *               std::vector<Tensor<2,dim> > &grads,
+ *               std::vector<Tensor<3,dim> > &grad_grads) const;
+ * @endcode
+ *
+ * The polynomial spaces are usually described as direct sum of
+ * simpler spaces. In such a case, the usual basis of node functionals
+ * is not dual to the basis of the polynomial space. Therefore, the
+ * matrix #inverse_node_matrix can be filled by the constructor of a
+ * derived class such that the usual interpolation condition
+ * <i>N<sub>i</sub>(v<sub>j</sub>)</i> holds on the reference cell.
+ *
+ * In many cases, the node functionals depend on the shape of the mesh
+ * cell, since they evaluate normal or tangential components on the
+ * faces. In order to allow for a set of transformations, the variable
+ * #mapping_type has been introduced. It should also be set in the
+ * constructor of a derived class.
+ *
+ * This class is not a fully implemented FiniteElement class, but
+ * implements some common features of vector valued elements based on
+ * vector valued polynomial classes. What's missing here in particular
+ * is information on the topological location of the node values.
+ *
+ * @see PolynomialsBDM, PolynomialsRaviartThomas
+ *
+ * @author Guido Kanschat, 2005
+ **/
+template <class POLY, int dim>
+class FE_PolyTensor : public FiniteElement<dim>
+{
+  public:
+                                    /**
+                                     * Constructor.
+                                     */
+    FE_PolyTensor (unsigned int degree,
+                  const FiniteElementData<dim> &fe_data,
+                  const std::vector<bool> &restriction_is_additive_flags,
+                  const std::vector<std::vector<bool> > &nonzero_components);
+
+                                    /**
+                                     * Since these elements are
+                                     * vector valued, an exception is
+                                     * thrown.
+                                     */
+    virtual double shape_value (const unsigned int i,
+                               const Point<dim> &p) const;
+    
+    virtual double shape_value_component (const unsigned int i,
+                                         const Point<dim> &p,
+                                         const unsigned int component) const;
+
+                                    /**
+                                     * Since these elements are
+                                     * vector valued, an exception is
+                                     * thrown.
+                                     */
+    virtual Tensor<1,dim> shape_grad (const unsigned int  i,
+                                     const Point<dim>   &p) const;
+
+    virtual Tensor<1,dim> shape_grad_component (const unsigned int i,
+                                               const Point<dim> &p,
+                                               const unsigned int component) const;
+
+                                    /**
+                                     * Since these elements are
+                                     * vector valued, an exception is
+                                     * thrown.
+                                     */
+    virtual Tensor<2,dim> shape_grad_grad (const unsigned int  i,
+                                          const Point<dim> &p) const;
+
+    virtual Tensor<2,dim> shape_grad_grad_component (const unsigned int i,
+                                                    const Point<dim> &p,
+                                                    const unsigned int component) const;
+
+                                     /**
+                                     * Number of base elements in a
+                                     * mixed discretization. Since
+                                     * this is not a composed element,
+                                     * return one.
+                                     */
+    virtual unsigned int n_base_elements () const;
+    
+                                    /**
+                                     * Access to base element
+                                     * objects. Since this element is
+                                     * not composed of several elements,
+                                     * <tt>base_element(0)</tt> is
+                                     * <tt>this</tt>, and all other
+                                     * indices throw an error.
+                                     */
+    virtual const FiniteElement<dim> &
+    base_element (const unsigned int index) const;
+
+                                     /**
+                                      * Multiplicity of base element
+                                      * <tt>index</tt>. Since this is
+                                      * not a composed element,
+                                      * <tt>element_multiplicity(0)</tt>
+                                      * returns one, and all other
+                                      * indices will throw an error.
+                                      */
+    virtual unsigned int element_multiplicity (const unsigned int index) const;
+
+    
+  protected:
+                                    /**
+                                     * Different options for
+                                     * transforming the basis
+                                     * functions from the reference
+                                     * cell to the actual mesh cell.
+                                     *
+                                     * Most vector valued elements
+                                     * either transform shape
+                                     * functions to keep node values
+                                     * on edges meaningful. Still, in
+                                     * special cases, it may be
+                                     * possible to avoid the mapping.
+                                     */
+    enum MappingType {
+                                          /// Shape functions do not depend on actual mesh cell
+         independent,
+                                          /// Shape functions are transformed covariant.
+         covariant,
+                                          /// Shape functions are transformed contravariant.
+         contravariant
+    };
+
+                                    /**
+                                     * The mapping type to be used to
+                                     * map shape functions from the
+                                     * reference cell to the msh
+                                     * cell.
+                                     */
+    MappingType mapping_type;
+    
+    virtual
+    typename Mapping<dim>::InternalDataBase *
+    get_data (const UpdateFlags,
+             const Mapping<dim>& mapping,
+             const Quadrature<dim>& quadrature) const ;
+
+    virtual void
+    fill_fe_values (const Mapping<dim> &mapping,
+                   const typename Triangulation<dim>::cell_iterator &cell,
+                   const Quadrature<dim>                &quadrature,
+                   typename Mapping<dim>::InternalDataBase      &mapping_internal,
+                   typename Mapping<dim>::InternalDataBase      &fe_internal,
+                   FEValuesData<dim>& data) const;
+    
+    virtual void
+    fill_fe_face_values (const Mapping<dim> &mapping,
+                        const typename Triangulation<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 ;
+    
+    virtual void
+    fill_fe_subface_values (const Mapping<dim> &mapping,
+                           const typename Triangulation<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 ;
+
+                                    /**
+                                     * Fields of cell-independent
+                                     * data for FE_PolyTensor. Stores
+                                     * the values of the shape
+                                     * functions and their
+                                     * derivatives on the reference
+                                     * cell for later use.
+                                     *
+                                     * All tables are organized in a
+                                     * way, that the value for shape
+                                     * function <i>i</i> at
+                                     * quadrature point <i>k</i> is
+                                     * accessed by indices
+                                     * <i>(i,k)</i>.
+                                     */
+    class InternalData : public FiniteElementBase<dim>::InternalDataBase
+    {
+      public:
+                                        /**
+                                         * Array with shape function
+                                         * values in quadrature
+                                         * points. There is one
+                                         * row for each shape
+                                         * function, containing
+                                         * values for each quadrature
+                                         * point.
+                                         */
+       std::vector<std::vector<Tensor<1,dim> > > shape_values;
+
+                                        /**
+                                         * Array with shape function
+                                         * gradients in quadrature
+                                         * points. There is one
+                                         * row for each shape
+                                         * function, containing
+                                         * values for each quadrature
+                                         * point.
+                                         */      
+       std::vector<std::vector<Tensor<2,dim> > > shape_grads;
+    };
+    
+                                    /**
+                                     * Degree of the polynomials.
+                                     */  
+    unsigned int degree;
+
+                                     /**
+                                      * The polynomial space. Its type
+                                      * is given by the template
+                                      * parameter POLY.
+                                      */    
+    POLY poly_space;
+                                    /**
+                                     * The inverse of the matrix
+                                     * <i>a<sub>ij</sub></i> of node
+                                     * values <i>N<sub>i</sub></i>
+                                     * applied to polynomial
+                                     * <i>p<sub>j</sub></i>. This
+                                     * matrix is used to convert
+                                     * polynomials in the "raw" basis
+                                     * provided in #poly_space to the
+                                     * basis dual to the node
+                                     * functionals on the reference cell.
+                                     *
+                                     * This object is not filled by
+                                     * FE_PolyTensor, but is a chance
+                                     * for a derived class to allow
+                                     * for reorganization of the
+                                     * basis functions. If it is left
+                                     * empty, the basis in
+                                     * #poly_space is used.
+                                     */
+    FullMatrix<double> inverse_node_matrix;
+
+                                    /**
+                                     * If a shape function is
+                                     * computed at a single point, we
+                                     * must compute all of them to
+                                     * apply #inverse_node_matrix. In
+                                     * order to avoid too much
+                                     * overhead, we cache the point
+                                     * and the function values for
+                                     * the next evaluation.
+                                     */ 
+    mutable Point<dim> cached_point;
+    
+                                    /**
+                                     * Cached shape function values after
+                                     * call to
+                                     * shape_value_component().
+                                     */
+    mutable std::vector<Tensor<1,dim> > cached_values;
+    
+                                    /**
+                                     * Cached shape function gradients after
+                                     * call to
+                                     * shape_grad_component().
+                                     */
+    mutable std::vector<Tensor<2,dim> > cached_grads;
+    
+                                    /**
+                                     * Cached second derivatives of
+                                     * shape functions after call to
+                                     * shape_grad_grad_component().
+                                     */
+    mutable std::vector<Tensor<3,dim> > cached_grad_grads;
+};
+
+/*@}*/
+
+#endif
diff --git a/deal.II/deal.II/source/fe/fe_poly_tensor.cc b/deal.II/deal.II/source/fe/fe_poly_tensor.cc
new file mode 100644 (file)
index 0000000..4c9cebd
--- /dev/null
@@ -0,0 +1,454 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2005 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/polynomials_bdm.h>
+#include <base/polynomials_raviart_thomas.h>
+#include <fe/fe_poly_tensor.h>
+#include <fe/fe_values.h>
+
+
+template <class POLY, int dim>
+FE_PolyTensor<POLY,dim>::FE_PolyTensor (
+  unsigned int degree,
+  const FiniteElementData<dim> &fe_data,
+  const std::vector<bool> &restriction_is_additive_flags,
+  const std::vector<std::vector<bool> > &nonzero_components):
+               FiniteElement<dim> (fe_data,
+                                   restriction_is_additive_flags,
+                                   nonzero_components),
+                degree(degree),
+                poly_space(POLY(degree))
+{
+  cached_point(0) = -1;
+}
+
+
+template <class POLY, int dim>
+double
+FE_PolyTensor<POLY,dim>::shape_value (
+  const unsigned int, const Point<dim> &) const
+{
+  Assert(false, typename FiniteElementBase<dim>::ExcFENotPrimitive());
+  return 0.;
+}
+
+
+template <class POLY, int dim>
+double
+FE_PolyTensor<POLY,dim>::shape_value_component (
+  const unsigned int i,
+  const Point<dim> &p,
+  const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component < dim, ExcIndexRange (component, 0, dim));
+  
+  if (cached_point != p || cached_values.size() == 0)
+    {
+      cached_point = p;
+      cached_values.resize(poly_space.n());
+      poly_space.compute(p, cached_values, cached_grads, cached_grad_grads);
+    }
+  
+  double s = 0;
+  if (inverse_node_matrix.n_cols() == 0)
+    return cached_values[i][component];
+  else
+    for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
+      s += inverse_node_matrix(i,j) * cached_values[j][component];
+  return s;
+}
+
+
+
+template <class POLY, int dim>
+Tensor<1,dim>
+FE_PolyTensor<POLY,dim>::shape_grad (
+  const unsigned int, const Point<dim> &) const
+{
+  Assert(false, typename FiniteElementBase<dim>::ExcFENotPrimitive());
+  return Tensor<1,dim>();
+}
+
+
+
+template <class POLY, int dim>
+Tensor<1,dim>
+FE_PolyTensor<POLY,dim>::shape_grad_component (
+  const unsigned int i,
+  const Point<dim> &p,
+  const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component < dim, ExcIndexRange (component, 0, dim));
+  
+  if (cached_point != p || cached_grads.size() == 0)
+    {
+      cached_point = p;
+      cached_grads.resize(poly_space.n());
+      poly_space.compute(p, cached_values, cached_grads, cached_grad_grads);
+    }
+  
+  Tensor<1,dim> s;
+  if (inverse_node_matrix.n_cols() == 0)
+    return cached_grads[i][component];
+  else
+    for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
+      s += inverse_node_matrix(i,j) * cached_grads[j][component];
+  
+  return s;
+}
+
+
+
+template <class POLY, int dim>
+Tensor<2,dim>
+FE_PolyTensor<POLY,dim>::shape_grad_grad (
+  const unsigned int, const Point<dim> &) const
+{
+  Assert(false, typename FiniteElementBase<dim>::ExcFENotPrimitive());
+  return Tensor<2,dim>();
+}
+
+
+
+template <class POLY, int dim>
+Tensor<2,dim>
+FE_PolyTensor<POLY,dim>::shape_grad_grad_component (
+  const unsigned int i,
+  const Point<dim> &p,
+  const unsigned int component) const
+{
+  Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
+  Assert (component < dim, ExcIndexRange (component, 0, dim));
+  
+  if (cached_point != p || cached_grad_grads.size() == 0)
+    {
+      cached_point = p;
+      cached_grad_grads.resize(poly_space.n());
+      poly_space.compute(p, cached_values, cached_grads, cached_grad_grads);
+    }
+  
+  Tensor<2,dim> s;
+  if (inverse_node_matrix.n_cols() == 0)
+    return cached_grad_grads[i][component];
+  else
+    for (unsigned int j=0;j<inverse_node_matrix.n_cols();++j)
+      s += inverse_node_matrix(i,j) * cached_grad_grads[j][component];
+  
+  return s;
+}
+
+
+
+//---------------------------------------------------------------------------
+// Data field initialization
+//---------------------------------------------------------------------------
+
+template <class POLY, int dim>
+typename Mapping<dim>::InternalDataBase *
+FE_PolyTensor<POLY,dim>::get_data (const UpdateFlags      update_flags,
+                            const Mapping<dim>    &mapping,
+                            const Quadrature<dim> &quadrature) const
+{
+                                  // generate a new data object and
+                                  // initialize some fields
+  InternalData* data = new InternalData;
+
+                                  // 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;
+
+                                  // some scratch arrays
+  std::vector<Tensor<1,dim> > values(0);
+  std::vector<Tensor<2,dim> > grads(0);
+  std::vector<Tensor<3,dim> > grad_grads(0);
+
+                                  // initialize fields only if really
+                                  // necessary. otherwise, don't
+                                  // allocate memory
+  if (flags & update_values)
+    {
+      values.resize (this->dofs_per_cell);
+      data->shape_values.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+       data->shape_values[i].resize (n_q_points);
+    }
+
+  if (flags & update_gradients)
+    {
+      grads.resize (this->dofs_per_cell);
+      data->shape_grads.resize (this->dofs_per_cell);
+      for (unsigned int i=0;i<this->dofs_per_cell;++i)
+       data->shape_grads[i].resize (n_q_points);
+    }
+
+                                  // if second derivatives through
+                                  // finite differencing is required,
+                                  // then initialize some objects for
+                                  // that
+  if (flags & update_second_derivatives)
+    {
+//      grad_grads.resize (this->dofs_per_cell);      
+      data->initialize_2nd (this, mapping, quadrature);
+    }
+
+                                  // Compute shape function values
+                                  // and derivatives on the reference
+                                  // cell. Make sure, that for the
+                                  // node values N_i holds
+                                  // N_i(v_j)=\delta_ij for all basis
+                                  // functions v_j
+  if (flags & (update_values | update_gradients))
+    for (unsigned int k=0; k<n_q_points; ++k)
+      {
+       poly_space.compute(quadrature.point(k),
+                          values, grads, grad_grads);
+       
+       if (flags & update_values)
+         if (inverse_node_matrix.n_cols() == 0)
+           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+             data->shape_values[i][k] = values[i];
+         else
+           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+             for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+               data->shape_values[i][k] = inverse_node_matrix(i,j) * values[j];
+       
+       if (flags & update_gradients)
+         if (inverse_node_matrix.n_cols() == 0)
+           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+             data->shape_grads[i][k] = grads[i];
+         else
+           for (unsigned int i=0; i<this->dofs_per_cell; ++i)
+             for (unsigned int j=0; j<this->dofs_per_cell; ++j)
+               data->shape_grads[i][k] = inverse_node_matrix(i,j) * grads[j];
+      }
+  return data;
+}
+
+
+
+
+//---------------------------------------------------------------------------
+// Fill data of FEValues
+//---------------------------------------------------------------------------
+
+template <class POLY, int dim>
+void
+FE_PolyTensor<POLY,dim>::fill_fe_values (const Mapping<dim>                   &mapping,
+                                  const typename Triangulation<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 unsigned int n_dofs = this->dofs_per_cell;
+  const unsigned int n_quad = quadrature.n_quadrature_points;
+  const UpdateFlags flags(fe_data.current_update_flags());
+  
+  Assert(mapping_type == independent, ExcNotImplemented());
+  
+  Assert(!(flags & update_values) || fe_data.shape_values.size() == n_dofs,
+        ExcDimensionMismatch(fe_data.shape_values.size(), n_dofs));
+  Assert(!(flags & update_values) || fe_data.shape_values[0].size() == n_quad,
+        ExcDimensionMismatch(fe_data.shape_values[0].size(), n_quad));
+  
+  for (unsigned int i=0; i<n_dofs; ++i)
+    {
+      const unsigned int first = data.shape_function_to_row_table[i];
+      
+      if (flags & update_values)
+       for (unsigned int k=0; k<n_quad; ++k)
+         for (unsigned int d=0;d<dim;++d)
+         data.shape_values(first+d,k) = fe_data.shape_values[i][k][d];
+      
+      if (flags & update_gradients)
+       {
+         std::vector<Tensor<2,dim> > shape_grads1 (n_quad);
+         mapping.transform_covariant(fe_data.shape_grads[i], 0,
+                                     shape_grads1,
+                                     mapping_data);
+         for (unsigned int k=0; k<n_quad; ++k)
+           for (unsigned int d=0;d<dim;++d)
+             data.shape_gradients[first+d][k] = shape_grads1[k][d];
+       }
+    }
+
+  const typename QProjector<dim>::DataSetDescriptor dsd;
+  if (flags & update_second_derivatives)
+    this->compute_2nd (mapping, cell, dsd.cell(),
+                      mapping_data, fe_data, data);
+}
+
+
+
+template <class POLY, int dim>
+void
+FE_PolyTensor<POLY,dim>::fill_fe_face_values (const Mapping<dim>                   &mapping,
+                               const typename Triangulation<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);
+
+  const unsigned int n_dofs = this->dofs_per_cell;
+  const unsigned int n_quad = quadrature.n_quadrature_points;
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                  // faces are stored contiguously)
+  
+  const typename QProjector<dim>::DataSetDescriptor dsd;
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = dsd.face (face, cell->face_orientation(face),
+               n_quad);
+  
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  Assert(mapping_type == independent, ExcNotImplemented());
+//TODO: Size assertions
+  
+  for (unsigned int i=0; i<n_dofs; ++i)
+    {
+      const unsigned int first = data.shape_function_to_row_table[i];
+      
+      if (flags & update_values)
+        for (unsigned int k=0; k<n_quad; ++k)
+         for (unsigned int d=0;d<dim;++d)
+           data.shape_values(first+d,k) = fe_data.shape_values[i][k+offset][d];
+      
+      if (flags & update_gradients)
+       {
+         std::vector<Tensor<2,dim> > shape_grads1 (n_quad);
+         mapping.transform_covariant(fe_data.shape_grads[i], offset,
+                                     shape_grads1, mapping_data);
+         for (unsigned int k=0; k<n_quad; ++k)
+           for (unsigned int d=0;d<dim;++d)
+             data.shape_gradients[first+d][k] = shape_grads1[k][d];
+       }
+    }
+
+  if (flags & update_second_derivatives)
+    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+}
+
+
+
+template <class POLY, int dim>
+void
+FE_PolyTensor<POLY,dim>::fill_fe_subface_values (const Mapping<dim>                   &mapping,
+                                  const typename Triangulation<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);
+
+  const unsigned int n_dofs = this->dofs_per_cell;
+  const unsigned int n_quad = quadrature.n_quadrature_points;
+                                  // offset determines which data set
+                                  // to take (all data sets for all
+                                // sub-faces are stored contiguously)
+
+  const typename QProjector<dim>::DataSetDescriptor dsd;
+  const typename QProjector<dim>::DataSetDescriptor offset
+    = dsd.sub_face (face, subface, cell->face_orientation(face),
+                   n_quad);
+
+  const UpdateFlags flags(fe_data.update_once | fe_data.update_each);
+
+  Assert(mapping_type == independent, ExcNotImplemented());
+//TODO: Size assertions
+  
+  for (unsigned int i=0; i<n_dofs; ++i)
+    {
+      const unsigned int first = data.shape_function_to_row_table[i];
+      
+      if (flags & update_values)
+        for (unsigned int k=0; k<n_quad; ++k)
+         for (unsigned int d=0;d<dim;++d)
+           data.shape_values(first+d,k) = fe_data.shape_values[i][k+offset][d];
+      
+      if (flags & update_gradients)
+       {
+         std::vector<Tensor<2,dim> > shape_grads1 (n_quad);
+         mapping.transform_covariant(fe_data.shape_grads[i], offset,
+                                     shape_grads1, mapping_data);
+         for (unsigned int k=0; k<n_quad; ++k)
+           for (unsigned int d=0;d<dim;++d)
+             data.shape_gradients[first+d][k] = shape_grads1[k][d];
+       }
+    }
+  
+  if (flags & update_second_derivatives)
+    this->compute_2nd (mapping, cell, offset, mapping_data, fe_data, data);
+}
+
+
+
+template <class POLY, int dim>
+unsigned int
+FE_PolyTensor<POLY,dim>::n_base_elements () const
+{
+  return 1;
+}
+
+
+
+template <class POLY, int dim>
+const FiniteElement<dim> &
+FE_PolyTensor<POLY,dim>::base_element (const unsigned int index) const
+{
+  Assert (index==0, ExcIndexRange(index, 0, 1));
+  return *this;
+}
+
+
+
+template <class POLY, int dim>
+unsigned int
+FE_PolyTensor<POLY,dim>::element_multiplicity (const unsigned int index) const
+{
+  Assert (index==0, ExcIndexRange(index, 0, 1));
+  return 1;
+}
+
+
+
+template class FE_PolyTensor<PolynomialsRaviartThomas<deal_II_dimension>,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.