]> https://gitweb.dealii.org/ - dealii.git/commitdiff
new wonderful maping
authorMarco Tezzele <marcotez@gmail.com>
Thu, 2 Apr 2015 12:48:28 +0000 (14:48 +0200)
committerMarco Tezzele <marcotez@gmail.com>
Sun, 5 Apr 2015 15:51:31 +0000 (17:51 +0200)
include/deal.II/fe/mapping_fe.h [new file with mode: 0644]
source/fe/CMakeLists.txt
source/fe/mapping_fe.cc [new file with mode: 0644]
source/fe/mapping_fe.inst.in [new file with mode: 0644]
tests/fe/mapping_fe_real_to_unit_q1.cc [new file with mode: 0644]
tests/fe/mapping_fe_real_to_unit_q1.output [new file with mode: 0644]
tests/fe/mapping_fe_real_to_unit_q5_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_real_to_unit_q5_curved.output [new file with mode: 0644]

diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h
new file mode 100644 (file)
index 0000000..a5e21c1
--- /dev/null
@@ -0,0 +1,742 @@
+// ---------------------------------------------------------------------
+// $Id: mapping_fe.h 30450 2013-08-23 15:48:29Z kronbichler $
+//
+// Copyright (C) 2001 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__mapping_fe_h
+#define __deal2__mapping_fe_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/table.h>
+#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/thread_management.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+/*!@addtogroup mapping */
+/*@{*/
+
+/**
+ * The MappingFE is a generalization of the MappingQEulerian class, for arbitrary
+ * vectorial finite elements. The main difference is that this class uses a vector
+ * of absolute positions instead of a vector of displacement.
+ * In particular we think of a collections of a FE_Q or
+ * Bezier finite element (FE_Bernstein) repeated a number of times equal to the space
+ * dimension. The idea is to construct the mapping using a vector of control
+ * points, a DoFHandler associated to the geometry of the problem and a
+ * ComponentMask that tells us which components to use for the mapping.
+ * This mapping will grab from the DoFHandler the finite element, or better
+ * the collection of finite elements, to compute the mapping shape functions.
+ * So we will have two different Finite Element and DoFHandler, one for the
+ * solution field and one to describe the geometry of the problem. Historically
+ * in the deal.II library there was not such a distinction. The differences
+ * between this mapping and the MappingQ class are quite important.
+ * The MappingFE, being a generalization, requires a higher level of abstraction.
+ * This is the reason why it takes a DoFHandler and a vector of control points
+ * that are the coefficients of the shape function (so in general it is a vector
+ * of coefficient).
+ *
+ *
+ * Typically, the DoFHandler operates on a finite element that is constructed
+ * as a system element (FESystem) from continuous FE_Q() objects. An example
+ * is shown below:
+ * @code
+ *    const FE_Q<dim,spacedim> feq(1);
+ *    const FESystem<dim,spacedim> fesystem(feq, spacedim);
+ *    DoFHandler<dim,spacedim> dhq(triangulation);
+ *    dhq.distribute_dofs(fesystem);
+ *    Vector<double> eulerq(dhq.n_dofs());
+ *    const ComponentMask mask(spacedim, true);
+ *    MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+ *    map.update_euler_vector_using_triangulation(eulerq);
+ * @endcode
+
+
+ *
+ * @author Luca Heltai, Marco Tezzele 2013, 2015
+ */
+template <int dim, int spacedim=dim,
+          class DH=DoFHandler<dim,spacedim>,
+          class VECTOR=Vector<double> >
+class MappingFE : public Mapping<dim,spacedim>
+{
+public:
+  /**
+   * Constructor. The first argument is a VECTOR that specifies the
+   * transformation of the domain from the reference to the current
+   * configuration. This is filled calling the method
+   * update_euler_vector_using_triangulation.
+   */
+  MappingFE (const VECTOR  &euler_vector,
+             const DH      &euler_dof_handler,
+             const ComponentMask mask=ComponentMask());
+
+  /**
+   * Copy constructor. Performs a deep copy, i.e. duplicates what #tensor_pols
+   * points to instead of simply copying the #tensor_pols pointer as done by a
+   * default copy constructor.
+   */
+  MappingFE (const MappingFE<dim,spacedim,DH,VECTOR> &mapping);
+
+  /**
+   * Destructor.
+   */
+  virtual ~MappingFE ();
+
+  /** Fill the euler vector with
+  the information coming from
+  the triangulation. Makes this
+  map equivalent to MappingQ1,
+  and it works ONLY if the
+  underlying fe has support
+  points. */
+  void update_euler_vector_using_triangulation(VECTOR &vector);
+
+
+
+  /**
+   * Transforms the point @p p on the unit cell to the point @p p_real on the
+   * real cell @p cell and returns @p p_real.
+   */
+  virtual Point<spacedim>
+  transform_unit_to_real_cell (
+    const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+    const Point<dim>                                 &p) const;
+
+  /**
+   * Transforms the point @p p on the real cell to the point @p p_unit on the
+   * unit cell @p cell and returns @p p_unit.
+   *
+   * Uses Newton iteration and the @p transform_unit_to_real_cell function.
+   *
+   * In the codimension one case, this function returns the normal projection
+   * of the real point @p p on the curve or surface identified by the @p cell.
+   */
+  virtual Point<dim>
+  transform_real_to_unit_cell (
+    const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+    const Point<spacedim>                            &p) const;
+
+
+  virtual void
+  transform (const VectorSlice<const std::vector<Tensor<1,dim> > > input,
+             VectorSlice<std::vector<Tensor<1,spacedim> > > output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+  virtual void
+  transform (const VectorSlice<const std::vector<DerivativeForm<1, dim, spacedim> > >    input,
+             VectorSlice<std::vector<Tensor<2,spacedim> > > output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+  virtual
+  void
+  transform (const VectorSlice<const std::vector<Tensor<2, dim> > >     input,
+             VectorSlice<std::vector<Tensor<2,spacedim> > >             output,
+             const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+             const MappingType type) const;
+
+
+
+  /**
+   * Return the degree of the mapping, i.e. the value which was passed to the
+   * constructor.
+   */
+  unsigned int get_degree () const;
+
+  /**
+   * Return the ComponentMask of the mapping, i.e. which components to use for
+   * the mapping.
+   */
+  ComponentMask get_fe_mask () const;
+
+  /**
+   * Return a pointer to a copy of the present object. The caller of this copy
+   * then assumes ownership of it.
+   */
+  virtual
+  Mapping<dim,spacedim> *clone () const;
+
+
+  /**
+   * Storage for internal data of
+   * d-linear transformation.
+   */
+  class InternalData : public Mapping<dim,spacedim>::InternalDataBase
+  {
+  public:
+    /**
+     * Constructor.
+     */
+    InternalData(const FiniteElement<dim,spacedim> &fe,
+                 const ComponentMask mask);
+
+    /**
+     * Shape function at quadrature
+     * point. Shape functions are
+     * in tensor product order, so
+     * vertices must be reordered
+     * to obtain transformation.
+     */
+    double shape (const unsigned int qpoint,
+                  const unsigned int shape_nr) const;
+
+    /**
+     * Shape function at quadrature
+     * point. See above.
+     */
+    double &shape (const unsigned int qpoint,
+                   const unsigned int shape_nr);
+
+    /**
+     * Gradient of shape function
+     * in quadrature point. See
+     * above.
+     */
+    Tensor<1,dim> derivative (const unsigned int qpoint,
+                              const unsigned int shape_nr) const;
+
+    /**
+     * Gradient of shape function
+     * in quadrature point. See
+     * above.
+     */
+    Tensor<1,dim> &derivative (const unsigned int qpoint,
+                               const unsigned int shape_nr);
+
+    /**
+     * Second derivative of shape
+     * function in quadrature
+     * point. See above.
+     */
+    Tensor<2,dim> second_derivative (const unsigned int qpoint,
+                                     const unsigned int shape_nr) const;
+
+    /**
+     * Second derivative of shape
+     * function in quadrature
+     * point. See above.
+     */
+    Tensor<2,dim> &second_derivative (const unsigned int qpoint,
+                                      const unsigned int shape_nr);
+
+    /**
+     * Return an estimate (in
+     * bytes) or the memory
+     * consumption of this
+     * object.
+     */
+    virtual std::size_t memory_consumption () const;
+
+    /**
+     * Values of shape
+     * functions. Access by
+     * function @p shape.
+     *
+     * Computed once.
+     */
+    std::vector<double> shape_values;
+
+    /**
+     * Values of shape function
+     * derivatives. Access by
+     * function @p derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<1,dim> > shape_derivatives;
+
+    /**
+     * Values of shape function
+     * second derivatives. Access
+     * by function
+     * @p second_derivative.
+     *
+     * Computed once.
+     */
+    std::vector<Tensor<2,dim> > shape_second_derivatives;
+
+    /**
+     * Tensors of covariant
+     * transformation at each of
+     * the quadrature points. The
+     * matrix stored is the
+     * Jacobian * G^{-1},
+     * where G = Jacobian^{t} * Jacobian,
+     * is the first fundamental
+     * form of the map;
+     * if dim=spacedim then
+     * it reduces to the transpose of the
+     * inverse of the Jacobian
+     * matrix, which itself is
+     * stored in the
+     * @p contravariant field of
+     * this structure.
+     *
+     * Computed on each cell.
+     */
+    std::vector<DerivativeForm<1,dim, spacedim > >  covariant;
+
+    /**
+     * Tensors of contravariant
+     * transformation at each of
+     * the quadrature points. The
+     * contravariant matrix is
+     * the Jacobian of the
+     * transformation,
+     * i.e. $J_{ij}=dx_i/d\hat x_j$.
+     *
+     * Computed on each cell.
+     */
+    std::vector< DerivativeForm<1,dim,spacedim> > contravariant;
+
+    /**
+     * Unit tangential vectors. Used
+     * for the computation of
+     * boundary forms and normal
+     * vectors.
+     *
+     * This vector has
+     * (dim-1)GeometryInfo::faces_per_cell
+     * entries. The first
+     * GeometryInfo::faces_per_cell
+     * contain the vectors in the first
+     * tangential direction for each
+     * face; the second set of
+     * GeometryInfo::faces_per_cell
+     * entries contain the vectors in the
+     * second tangential direction (only
+     * in 3d, since there we have 2
+     * tangential directions per face),
+     * etc.
+     *
+     * Filled once.
+     */
+    std::vector<std::vector<Tensor<1,dim> > > unit_tangentials;
+
+    /**
+     * Auxiliary vectors for internal use.
+     */
+    std::vector<std::vector<Tensor<1,spacedim> > > aux;
+
+    /**
+     * Number of shape
+     * functions. If this is a Q1
+     * mapping, then it is simply
+     * the number of vertices per
+     * cell. However, since also
+     * derived classes use this
+     * class (e.g. the
+     * Mapping_Q() class),
+     * the number of shape
+     * functions may also be
+     * different.
+     */
+    unsigned int n_shape_functions;
+
+    ComponentMask mask;
+  };
+
+
+  /**
+   * Transforms a point @p p on
+   * the unit cell to the point
+   * @p p_real on the real cell
+   * @p cell and returns @p p_real.
+   *
+   * This function is called by
+   * @p transform_unit_to_real_cell
+   * and multiple times (through the
+   * Newton iteration) by
+   * @p transform_real_to_unit_cell_internal.
+   *
+   * Takes a reference to an
+   * @p InternalData that must
+   * already include the shape
+   * values at point @p p and the
+   * mapping support points of the
+   * cell.
+   *
+   * This @p InternalData argument
+   * avoids multiple computations
+   * of the shape values at point
+   * @p p and especially multiple
+   * computations of the mapping
+   * support points.
+   */
+  Point<spacedim>
+  transform_unit_to_real_cell_internal (const InternalData &mdata) const;
+
+
+  /**
+   * Transforms the point @p p on
+   * the real cell to the corresponding
+   * point on the unit cell
+   * @p cell by a Newton
+   * iteration.
+   *
+   * Takes a reference to an
+   * @p InternalData that is
+   * assumed to be previously
+   * created by the @p get_data
+   * function with @p UpdateFlags
+   * including
+   * @p update_transformation_values
+   * and
+   * @p update_transformation_gradients
+   * and a one point Quadrature
+   * that includes the given
+   * initial guess for the
+   * transformation
+   * @p initial_p_unit.  Hence this
+   * function assumes that
+   * @p mdata already includes the
+   * transformation shape values
+   * and gradients computed at
+   * @p initial_p_unit.
+   *
+   * @p mdata will be changed by
+   * this function.
+   */
+  Point<dim>
+  transform_real_to_unit_cell_internal (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                        const Point<spacedim> &p,
+                                        const Point<dim> &initial_p_unit,
+                                        InternalData &mdata) const;
+
+  /**
+   * Do the computation for the
+   * <tt>fill_*</tt> functions.
+   */
+  void compute_fill (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                     const unsigned int      npts,
+                     const typename QProjector<dim>::DataSetDescriptor data_set,
+                     const CellSimilarity::Similarity cell_similarity,
+                     InternalData           &data,
+                     std::vector<Point<spacedim> > &quadrature_points) const;
+
+
+  /**
+   * Do the computation for the
+   * <tt>fill_*</tt> functions.
+   */
+  void compute_fill_face (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                          const unsigned int      face_no,
+                          const unsigned int      subface_no,
+                          const unsigned int      npts,
+                          const typename QProjector<dim>::DataSetDescriptor data_set,
+                          const std::vector<double>   &weights,
+                          InternalData           &mapping_data,
+                          std::vector<Point<spacedim> >    &quadrature_points,
+                          std::vector<double>         &JxW_values,
+                          std::vector<Tensor<1,spacedim> > &boundary_form,
+                          std::vector<Point<spacedim> > &normal_vectors,
+                          std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+                          std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const;
+
+
+  /**
+   * Always returns @p false.
+   */
+  virtual
+  bool preserves_vertex_locations () const;
+
+  DeclException0(ExcInactiveCell);
+
+protected:
+  /**
+   * Implementation of the interface in Mapping.
+   */
+  virtual void
+  fill_fe_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                  const Quadrature<dim>                                     &quadrature,
+                  typename Mapping<dim,spacedim>::InternalDataBase          &mapping_data,
+                  typename std::vector<Point<spacedim> >                    &quadrature_points,
+                  std::vector<double>                                       &JxW_values,
+                  std::vector<DerivativeForm<1,dim,spacedim> >       &jacobians,
+                  std::vector<DerivativeForm<2,dim,spacedim> >       &jacobian_grads,
+                  std::vector<DerivativeForm<1,spacedim,dim> >      &inverse_jacobians,
+                  std::vector<Point<spacedim> >                             &cell_normal_vectors,
+                  CellSimilarity::Similarity                           &cell_similarity) const ;
+
+  /**
+   * Implementation of the interface in Mapping.
+   */
+  virtual void
+  fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                       const unsigned int face_no,
+                       const Quadrature<dim-1>& quadrature,
+                       typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                       typename std::vector<Point<spacedim> >        &quadrature_points,
+                       std::vector<double>             &JxW_values,
+                       std::vector<Tensor<1,spacedim> >             &exterior_forms,
+                       std::vector<Point<spacedim> >                &normal_vectors,
+                       std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+                       std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const ;
+
+  /**
+   * Implementation of the interface in Mapping.
+   */
+  virtual void
+  fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                          const unsigned int face_no,
+                          const unsigned int sub_no,
+                          const Quadrature<dim-1>& quadrature,
+                          typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                          typename std::vector<Point<spacedim> >        &quadrature_points,
+                          std::vector<double>             &JxW_values,
+                          std::vector<Tensor<1,spacedim> > &exterior_forms,
+                          std::vector<Point<spacedim> >    &normal_vectors,
+                          std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+                          std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const ;
+
+
+  /**
+     This function and the next allow to generate the transform require by
+     the virtual transform() in mapping, but unfortunately in C++ one cannot
+     declare a virtual template function.
+  */
+  template < int rank >
+  void
+  transform_fields(const VectorSlice<const std::vector<Tensor<rank,dim>      > > input,
+                   VectorSlice<      std::vector<Tensor<rank,spacedim> > > output,
+                   const typename Mapping<dim,spacedim>::InternalDataBase &internal,
+                   const MappingType type) const;
+
+
+  /**
+       see doc in transform_fields
+   */
+  template < int rank >
+  void
+  transform_differential_forms(
+    const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > >    input,
+    VectorSlice<std::vector<DerivativeForm<rank, spacedim,spacedim> > > output,
+    const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+    const MappingType mapping_type) const;
+
+
+protected:
+  /**
+  * Reference to the vector of shifts.
+  */
+
+  SmartPointer<const VECTOR, MappingFE<dim,spacedim,DH,VECTOR> >euler_vector;
+  /**
+   * A FiniteElement object which is only needed in 3D, since it knows how to reorder
+   * shape functions/DoFs on non-standard faces. This is used to reorder
+   * support points in the same way. We could make this a pointer to prevent
+   * construction in 1D and 2D, but since memory and time requirements are not
+   * particularly high this seems unnecessary at the moment.
+   */
+  SmartPointer<const FiniteElement<dim,spacedim>, MappingFE<dim,spacedim,DH,VECTOR> > fe;
+
+
+  /**
+   * Pointer to the DoFHandler to which the mapping vector is associated.
+   */
+  SmartPointer<const DH,MappingFE<dim,spacedim,DH,VECTOR> >euler_dof_handler;
+
+
+
+private:
+//
+  /**
+   * Update internal degrees of
+   * freedom. */
+  void update_internal_dofs(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
+
+  mutable std::vector<double> local_dofs;
+
+  mutable std::vector<unsigned int> dof_indices;
+
+  /**
+   * Mutex to protect local_dofs.
+   */
+
+  mutable Threads::Mutex mutex;
+
+
+  virtual void
+  compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
+                          typename MappingFE<dim, spacedim>::InternalData &data) const;
+
+  UpdateFlags
+  update_once (const UpdateFlags in) const;
+
+  UpdateFlags
+  update_each (const UpdateFlags in) const;
+
+  void
+  compute_data (const UpdateFlags      update_flags,
+                const Quadrature<dim>  &q,
+                const unsigned int     n_original_q_points,
+                InternalData           &data) const;
+
+  void
+  compute_face_data (const UpdateFlags      update_flags,
+                     const Quadrature<dim>  &q,
+                     const unsigned int     n_original_q_points,
+                     InternalData           &data) const;
+
+  virtual
+  typename Mapping<dim,spacedim>::InternalDataBase *
+  get_data (const UpdateFlags,
+            const Quadrature<dim> &quadrature) const;
+
+  virtual
+  typename Mapping<dim,spacedim>::InternalDataBase *
+  get_face_data (const UpdateFlags flags,
+                 const Quadrature<dim-1>& quadrature) const;
+
+  virtual
+  typename Mapping<dim,spacedim>::InternalDataBase *
+  get_subface_data (const UpdateFlags flags,
+                    const Quadrature<dim-1>& quadrature) const;
+
+
+  /*
+   * Which components to use for the mapping.
+   */
+  const ComponentMask fe_mask;
+
+
+  /**
+   * Mapping between indices in the FE space and the real space. This vector contains one
+   * index for each component of the finite element space. If the index is one for which
+   * the ComponentMask which is used to construct this element is false, then
+   * numbers::invalid_unsigned_int is returned, otherwise the component in real space is
+   * returned. For example, if we construct the mapping using ComponentMask(spacedim, true),
+   * then this vector contains {0,1,2} in spacedim = 3.
+   */
+  std::vector<unsigned int> fe_to_real;
+
+
+  /**
+   * Declare other MappingFE classes friends.
+   */
+  template <int,int,class,class> friend class MappingFE;
+};
+
+/*@}*/
+
+/* -------------- declaration of explicit specializations ------------- */
+
+#ifndef DOXYGEN
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+inline
+double
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::shape (const unsigned int qpoint,
+                                                        const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_values.size()));
+  return shape_values [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+inline
+double &
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::shape (const unsigned int qpoint,
+                                                        const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_values.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_values.size()));
+  return shape_values [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+inline
+Tensor<1,dim>
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int qpoint,
+    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_derivatives.size()));
+  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+inline
+Tensor<1,dim> &
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int qpoint,
+    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_derivatives.size()));
+  return shape_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, class DH, class VECTOR>
+inline
+Tensor<2,dim>
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::second_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr) const
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_second_derivatives.size()));
+  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+
+template <int dim, int spacedim, class DH, class VECTOR>
+inline
+Tensor<2,dim> &
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::second_derivative (const unsigned int qpoint,
+    const unsigned int shape_nr)
+{
+  Assert(qpoint*n_shape_functions + shape_nr < shape_second_derivatives.size(),
+         ExcIndexRange(qpoint*n_shape_functions + shape_nr, 0,
+                       shape_second_derivatives.size()));
+  return shape_second_derivatives [qpoint*n_shape_functions + shape_nr];
+}
+
+
+template <int dim, int spacedim, class DH, class VECTOR>
+inline
+bool
+MappingFE<dim,spacedim,DH,VECTOR>::preserves_vertex_locations () const
+{
+  return false;
+}
+
+
+
+
+#endif // DOXYGEN
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index b676e48baf5f0dc21b9cf3bf13df6a25a9a69ea1..80cfbcaa600ae6cfcc27ce047fd9461f0dbf6e34 100644 (file)
@@ -48,6 +48,7 @@ SET(_src
   mapping_c1.cc
   mapping_cartesian.cc
   mapping.cc
+  mapping_fe.cc
   mapping_q1.cc
   mapping_q1_eulerian.cc
   mapping_q.cc
@@ -87,6 +88,7 @@ SET(_inst
   mapping_c1.inst.in
   mapping_cartesian.inst.in
   mapping.inst.in
+  mapping_fe.inst.in
   mapping_q1_eulerian.inst.in
   mapping_q1.inst.in
   mapping_q_eulerian.inst.in
diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc
new file mode 100644 (file)
index 0000000..e66dd6e
--- /dev/null
@@ -0,0 +1,1386 @@
+// ---------------------------------------------------------------------
+// $Id: mapping_fe.cc 30450 2013-08-23 15:48:29Z kronbichler $
+//
+// Copyright (C) 2001 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/polynomial.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/base/qprojector.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <numeric>
+#include <memory>
+#include <fstream>
+
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::InternalData (const FiniteElement<dim,spacedim> &fe,
+    const ComponentMask mask)
+  :
+  n_shape_functions (fe.dofs_per_cell),
+  mask (mask)
+{}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+std::size_t
+MappingFE<dim,spacedim,DH,VECTOR>::InternalData::memory_consumption () const
+{
+  return 0;
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const VECTOR  &euler_vector,
+                                              const DH      &euler_dof_handler,
+                                              const ComponentMask mask)
+  :
+  euler_vector(&euler_vector),
+  fe(&euler_dof_handler.get_fe()),
+  euler_dof_handler(&euler_dof_handler),
+  local_dofs(fe->dofs_per_cell),
+  dof_indices(fe->dofs_per_cell),
+  fe_mask(mask.size() ? mask :
+          ComponentMask(fe->get_nonzero_components(0).size(), true)),
+  fe_to_real(fe_mask.size(), numbers::invalid_unsigned_int)
+{
+
+  unsigned int size = 0;
+  for (unsigned int i=0; i<fe_mask.size(); ++i)
+    {
+      if (fe_mask[i])
+        fe_to_real[i] = size++;
+    }
+  AssertDimension(size,spacedim);
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const MappingFE<dim,spacedim,DH,VECTOR> &mapping)
+  :
+  euler_vector(mapping.euler_vector),
+  fe(mapping.fe),
+  euler_dof_handler(mapping.euler_dof_handler),
+  local_dofs(mapping.local_dofs),
+  dof_indices(mapping.dof_indices),
+  fe_mask(mapping.fe_mask),
+  fe_to_real(mapping.fe_to_real)
+{}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+MappingFE<dim,spacedim,DH,VECTOR>::~MappingFE ()
+{
+  euler_dof_handler = NULL;
+  fe = NULL;
+  euler_vector = NULL;
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::compute_shapes_virtual (
+  const std::vector<Point<dim> > &unit_points,
+  typename MappingFE<dim, spacedim>::InternalData &data) const
+{
+  const unsigned int n_points=unit_points.size();
+
+  if (data.shape_values.size()!=0 || data.shape_derivatives.size()!=0)
+    for (unsigned int point=0; point<n_points; ++point)
+      {
+        if (data.shape_values.size()!=0)
+          for (unsigned int i=0; i<data.n_shape_functions; ++i)
+            data.shape(point, i) = fe->shape_value(i, unit_points[point]);
+
+        if (data.shape_derivatives.size()!=0)
+          for (unsigned int i=0; i<data.n_shape_functions; ++i)
+            data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
+
+        if (data.shape_second_derivatives.size()!=0)
+          for (unsigned int i=0; i<data.n_shape_functions; ++i)
+            data.second_derivative(point, i) = fe->shape_grad_grad(i, unit_points[point]);
+      }
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+UpdateFlags
+MappingFE<dim,spacedim,DH,VECTOR>::update_once (const UpdateFlags in) const
+{
+  UpdateFlags out = UpdateFlags(in & (update_transformation_values
+                                      | update_transformation_gradients));
+
+  // Shape function values
+  if (in & update_quadrature_points)
+    out |= update_transformation_values;
+
+  // Shape function gradients
+  if (in & (update_covariant_transformation
+            | update_contravariant_transformation
+            | update_JxW_values
+            | update_boundary_forms
+            | update_normal_vectors
+            | update_jacobians
+            | update_jacobian_grads
+            | update_inverse_jacobians))
+    out |= update_transformation_gradients;
+
+  return out;
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+UpdateFlags
+MappingFE<dim,spacedim,DH,VECTOR>::update_each (const UpdateFlags in) const
+{
+  // Select flags of concern for the
+  // transformation.
+  UpdateFlags out = UpdateFlags(in & (update_quadrature_points
+                                      | update_covariant_transformation
+                                      | update_contravariant_transformation
+                                      | update_JxW_values
+                                      | update_boundary_forms
+                                      | update_normal_vectors
+                                      | update_volume_elements
+                                      | update_jacobians
+                                      | update_jacobian_grads
+                                      | update_inverse_jacobians));
+
+  // add flags if the respective
+  // quantities are necessary to
+  // compute what we need. note that
+  // some flags appear in both
+  // conditions and in subsequents
+  // set operations. this leads to
+  // some circular logic. the only
+  // way to treat this is to
+  // iterate. since there are 4
+  // if-clauses in the loop, it will
+  // take at most 3 iterations to
+  // converge. do them:
+  for (unsigned int i=0; i<4; ++i)
+    {
+      // The following is a little incorrect:
+      // If not applied on a face,
+      // update_boundary_forms does not
+      // make sense. On the other hand,
+      // it is necessary on a
+      // face. Currently,
+      // update_boundary_forms is simply
+      // ignored for the interior of a
+      // cell.
+      if (out & (update_JxW_values
+                 | update_normal_vectors))
+        out |= update_boundary_forms;
+
+      if (out & (update_covariant_transformation
+                 | update_JxW_values
+                 | update_jacobians
+                 | update_jacobian_grads
+                 | update_boundary_forms
+                 | update_normal_vectors))
+        out |= update_contravariant_transformation;
+
+      if (out & (update_inverse_jacobians))
+        out |= update_covariant_transformation;
+
+      // The contravariant transformation
+      // is a Piola transformation, which
+      // requires the determinant of the
+      // Jacobi matrix of the transformation.
+      // Therefore these values have to be
+      // updated for each cell.
+      if (out & update_contravariant_transformation)
+        out |= update_JxW_values;
+
+      if (out & update_normal_vectors)
+        out |= update_JxW_values;
+    }
+
+  return out;
+}
+
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::compute_data (const UpdateFlags      update_flags,
+                                                 const Quadrature<dim>  &q,
+                                                 const unsigned int     n_original_q_points,
+                                                 InternalData           &data) const
+{
+  const unsigned int n_q_points = q.size();
+
+  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);
+
+  if (flags & update_transformation_values)
+    data.shape_values.resize(data.n_shape_functions * n_q_points);
+
+  if (flags & update_transformation_gradients)
+    data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
+
+  if (flags & update_covariant_transformation)
+    data.covariant.resize(n_original_q_points);
+
+  if (flags & update_contravariant_transformation)
+    data.contravariant.resize(n_original_q_points);
+
+  if (flags & update_volume_elements)
+    data.volume_elements.resize(n_original_q_points);
+
+  if (flags & update_jacobian_grads)
+    data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
+
+  compute_shapes_virtual (q.get_points(), data);
+
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::compute_face_data (const UpdateFlags update_flags,
+                                                      const Quadrature<dim> &q,
+                                                      const unsigned int n_original_q_points,
+                                                      InternalData &data) const
+{
+  compute_data (update_flags, q, n_original_q_points, data);
+
+  if (dim > 1)
+    {
+      if (data.update_flags & update_boundary_forms)
+        {
+          data.aux.resize (dim-1, std::vector<Tensor<1,spacedim> > (n_original_q_points));
+
+          // Compute tangentials to the
+          // unit cell.
+          const unsigned int nfaces = GeometryInfo<dim>::faces_per_cell;
+          data.unit_tangentials.resize (nfaces*(dim-1),
+                                        std::vector<Tensor<1,dim> > (n_original_q_points));
+          if (dim==2)
+            {
+              // ensure a counterclock wise
+              // orientation of tangentials
+              static const int tangential_orientation[4]= {-1,1,1,-1};
+              for (unsigned int i=0; i<nfaces; ++i)
+                {
+                  Tensor<1,dim> tang;
+                  tang[1-i/2]=tangential_orientation[i];
+                  std::fill (data.unit_tangentials[i].begin(),
+                             data.unit_tangentials[i].end(), tang);
+                }
+            }
+          else if (dim==3)
+            {
+              for (unsigned int i=0; i<nfaces; ++i)
+                {
+                  Tensor<1,dim> tang1, tang2;
+
+                  const unsigned int nd=
+                    GeometryInfo<dim>::unit_normal_direction[i];
+
+                  // first tangential
+                  // vector in direction
+                  // of the (nd+1)%3 axis
+                  // and inverted in case
+                  // of unit inward normal
+                  tang1[(nd+1)%dim]=GeometryInfo<dim>::unit_normal_orientation[i];
+                  // second tangential
+                  // vector in direction
+                  // of the (nd+2)%3 axis
+                  tang2[(nd+2)%dim]=1.;
+
+                  // same unit tangents
+                  // for all quadrature
+                  // points on this face
+                  std::fill (data.unit_tangentials[i].begin(),
+                             data.unit_tangentials[i].end(), tang1);
+                  std::fill (data.unit_tangentials[nfaces+i].begin(),
+                             data.unit_tangentials[nfaces+i].end(), tang2);
+                }
+            }
+        }
+    }
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+typename Mapping<dim,spacedim>::InternalDataBase *
+MappingFE<dim,spacedim,DH,VECTOR>::get_data (const UpdateFlags update_flags,
+                                             const Quadrature<dim> &quadrature) const
+{
+  InternalData *data = new InternalData(*fe, fe_mask);
+  this->compute_data (update_flags, quadrature,
+                      quadrature.size(), *data);
+  return data;
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+typename Mapping<dim,spacedim>::InternalDataBase *
+MappingFE<dim,spacedim,DH,VECTOR>::get_face_data (const UpdateFlags update_flags,
+                                                  const Quadrature<dim-1>& quadrature) const
+{
+  InternalData *data = new InternalData(*fe, fe_mask);
+  const Quadrature<dim> q (QProjector<dim>::project_to_all_faces(quadrature));
+  this->compute_face_data (update_flags, q,
+                           quadrature.size(), *data);
+
+  return data;
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+typename Mapping<dim,spacedim>::InternalDataBase *
+MappingFE<dim,spacedim,DH,VECTOR>::get_subface_data (const UpdateFlags update_flags,
+                                                     const Quadrature<dim-1>& quadrature) const
+{
+  InternalData *data = new InternalData(*fe, fe_mask);
+  const Quadrature<dim> q (QProjector<dim>::project_to_all_subfaces(quadrature));
+  this->compute_face_data (update_flags, q,
+                           quadrature.size(), *data);
+
+  return data;
+}
+
+
+// Note that the CellSimilarity flag is modifyable, since MappingFE can need to
+// recalculate data even when cells are similar.
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_values (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+  const Quadrature<dim>                                     &q,
+  typename Mapping<dim,spacedim>::InternalDataBase          &mapping_data,
+  std::vector<Point<spacedim> >                             &quadrature_points,
+  std::vector<double>                                       &JxW_values,
+  std::vector<DerivativeForm<1,dim,spacedim>   >          &jacobians,
+  std::vector<DerivativeForm<2,dim,spacedim>    >           &jacobian_grads,
+  std::vector<DerivativeForm<1,spacedim,dim>    >           &inverse_jacobians,
+  std::vector<Point<spacedim> >                             &normal_vectors,
+  CellSimilarity::Similarity                                &cell_similarity) const
+{
+  AssertDimension(fe->dofs_per_cell, local_dofs.size());
+  Assert(local_dofs.size()>0, ExcMessage("Cannot do anything with zero degrees of freedom"));
+
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+  Assert (dynamic_cast<InternalData *> (&mapping_data) != 0, ExcInternalError());
+  InternalData &data = static_cast<InternalData &> (mapping_data);
+
+  // depending on this result, use this or the other data object for the
+  // mapping. furthermore, we need to ensure that the flag indicating whether
+  // we can use some similarity has to be modified - for a general MappingFE,
+  // the data needs to be recomputed anyway since then the mapping changes the
+  // data. this needs to be known also for later operations, so modify the
+  // variable here. this also affects the calculation of the next cell -- if
+  // we use Q1 data on the next cell, the data will still be invalid.
+
+  if (get_degree() > 1)
+    cell_similarity = CellSimilarity::invalid_next_cell;
+
+  const unsigned int n_q_points=q.size();
+
+  compute_fill (cell, n_q_points, QProjector<dim>::DataSetDescriptor::cell (), cell_similarity,
+                data, quadrature_points);
+
+  const UpdateFlags update_flags(data.current_update_flags());
+  const std::vector<double> &weights=q.get_weights();
+
+  // Multiply quadrature weights by absolute value of Jacobian determinants or
+  // the area element g=sqrt(DX^t DX) in case of codim > 0
+
+  if (update_flags & (update_normal_vectors | update_JxW_values))
+    {
+      AssertDimension (JxW_values.size(), n_q_points);
+
+      Assert( !(update_flags & update_normal_vectors ) ||
+              (normal_vectors.size() == n_q_points),
+              ExcDimensionMismatch(normal_vectors.size(), n_q_points));
+
+
+      if (cell_similarity != CellSimilarity::translation)
+        for (unsigned int point=0; point<n_q_points; ++point)
+          {
+            if (dim == spacedim)
+              {
+                const double det = data.contravariant[point].determinant();
+
+                // check for distorted cells.
+
+                // TODO: this allows for anisotropies of up to 1e6 in 3D and
+                // 1e12 in 2D. might want to find a finer
+                // (dimension-independent) criterion
+                Assert (det > 1e-12*Utilities::fixed_power<dim>(cell->diameter()/
+                                                                std::sqrt(double(dim))),
+                        (typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), det, point)));
+                JxW_values[point] = weights[point] * det;
+              }
+            // if dim==spacedim, then there is no cell normal to
+            // compute. since this is for FEValues (and not FEFaceValues),
+            // there are also no face normals to compute
+            else //codim>0 case
+              {
+                Tensor<1, spacedim> DX_t [dim];
+                for (unsigned int i=0; i<spacedim; ++i)
+                  for (unsigned int j=0; j<dim; ++j)
+                    DX_t[j][i] = data.contravariant[point][i][j];
+
+                Tensor<2, dim> G; //First fundamental form
+                for (unsigned int i=0; i<dim; ++i)
+                  for (unsigned int j=0; j<dim; ++j)
+                    G[i][j] = DX_t[i] * DX_t[j];
+
+                JxW_values[point] = sqrt(determinant(G)) * weights[point];
+
+                if (cell_similarity == CellSimilarity::inverted_translation)
+                  {
+                    // we only need to flip the normal
+                    if (update_flags & update_normal_vectors)
+                      normal_vectors[point] *= -1.;
+                  }
+                else
+                  {
+                    const unsigned int codim = spacedim-dim;
+
+                    if (update_flags & update_normal_vectors)
+                      {
+                        Assert( codim==1 , ExcMessage("There is no cell normal in codim 2."));
+
+                        if (dim==1)
+                          cross_product(normal_vectors[point], -DX_t[0]);
+                        else //dim == 2
+                          cross_product(normal_vectors[point],DX_t[0],DX_t[1]);
+
+                        normal_vectors[point] /= normal_vectors[point].norm();
+
+                        if (cell->direction_flag() == false)
+                          normal_vectors[point] *= -1.;
+                      }
+
+                  }
+              } //codim>0 case
+
+          }
+    }
+
+
+  // copy values from InternalData to vector given by reference
+  if (update_flags & update_jacobians)
+    {
+      AssertDimension (jacobians.size(), n_q_points);
+      if (cell_similarity != CellSimilarity::translation)
+        for (unsigned int point=0; point<n_q_points; ++point)
+          jacobians[point] = data.contravariant[point];
+    }
+
+
+  // calculate values of the derivatives of the Jacobians. do it here, since
+  // we only do it for cells, not faces.
+  if (update_flags & update_jacobian_grads)
+    {
+      AssertDimension (jacobian_grads.size(), n_q_points);
+
+      if (cell_similarity != CellSimilarity::translation)
+        {
+          std::fill(jacobian_grads.begin(),
+                    jacobian_grads.end(),
+                    DerivativeForm<2,dim,spacedim>());
+
+          const unsigned int data_set = QProjector<dim>::DataSetDescriptor::cell();
+
+          for (unsigned int point=0; point<n_q_points; ++point)
+            {
+              const Tensor<2,dim> *second =
+                &data.second_derivative(point+data_set, 0);
+
+              double result [spacedim][dim][dim];
+
+              for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                {
+                  unsigned int comp_k = fe->system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    for (unsigned int j=0; j<dim; ++j)
+                      for (unsigned int l=0; l<dim; ++l)
+                        result[fe_to_real[comp_k]][j][l] += (second[k][j][l]
+                                                             * local_dofs[k]);
+                }
+
+              // never touch any data for j=dim in case dim<spacedim, so it
+              // will always be zero as it was initialized
+              for (unsigned int i=0; i<spacedim; ++i)
+                for (unsigned int j=0; j<dim; ++j)
+                  for (unsigned int l=0; l<dim; ++l)
+                    jacobian_grads[point][i][j][l] = result[i][j][l];
+            }
+        }
+    }
+
+
+  // copy values from InternalData to vector given by reference
+  if (update_flags & update_inverse_jacobians)
+    {
+      AssertDimension (inverse_jacobians.size(), n_q_points);
+      if (cell_similarity != CellSimilarity::translation)
+        for (unsigned int point=0; point<n_q_points; ++point)
+          inverse_jacobians[point] = data.covariant[point].transpose();
+    }
+
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_face_values (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+  const unsigned int       face_no,
+  const Quadrature<dim-1> &q,
+  typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  std::vector<Point<spacedim> >     &quadrature_points,
+  std::vector<double>          &JxW_values,
+  std::vector<Tensor<1,spacedim> >             &exterior_forms,
+  std::vector<Point<spacedim> >                &normal_vectors,
+  std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+  std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const
+// std::vector<Tensor<1,spacedim> >  &exterior_forms,
+// std::vector<Point<spacedim> >     &normal_vectors) const
+{
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+
+//  AssertThrow(false, ExcNotImplemented());
+
+
+
+  Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
+          ExcInternalError());
+  InternalData &data = static_cast<InternalData &> (mapping_data);
+
+  const unsigned int n_q_points=q.size();
+  this->compute_fill_face (cell, face_no, numbers::invalid_unsigned_int,
+                           n_q_points,
+                           QProjector<dim>::DataSetDescriptor::
+                           face (face_no,
+                                 cell->face_orientation(face_no),
+                                 cell->face_flip(face_no),
+                                 cell->face_rotation(face_no),
+                                 n_q_points),
+                           q.get_weights(),
+                           data,
+                           quadrature_points, JxW_values,
+                           exterior_forms, normal_vectors, jacobians,
+                           inverse_jacobians);
+  // quadrature_points, JxW_values,
+  // exterior_forms, normal_vectors);
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_subface_values (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                                                           const unsigned int       face_no,
+                                                           const unsigned int       sub_no,
+                                                           const Quadrature<dim-1> &q,
+                                                           typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+                                                           std::vector<Point<spacedim> >     &quadrature_points,
+                                                           std::vector<double>          &JxW_values,
+                                                           std::vector<Tensor<1,spacedim> > &exterior_forms,
+                                                           std::vector<Point<spacedim> >    &normal_vectors,
+                                                           std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+                                                           std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const
+// std::vector<Tensor<1,spacedim> >  &exterior_forms,
+// std::vector<Point<spacedim> >     &normal_vectors) const
+{
+  //AssertThrow(false, ExcNotImplemented());
+
+
+  // convert data object to internal data for this class. fails with an
+  // exception if that is not possible
+  Assert (dynamic_cast<InternalData *> (&mapping_data) != 0,
+          ExcInternalError());
+  InternalData &data = static_cast<InternalData &> (mapping_data);
+
+  const unsigned int n_q_points=q.size();
+  this->compute_fill_face (cell, face_no, sub_no,
+                           n_q_points,
+                           QProjector<dim>::DataSetDescriptor::
+                           subface (face_no, sub_no,
+                                    cell->face_orientation(face_no),
+                                    cell->face_flip(face_no),
+                                    cell->face_rotation(face_no),
+                                    n_q_points,
+                                    cell->subface_case(face_no)),
+                           q.get_weights(),
+                           data,
+                           quadrature_points, JxW_values,
+                           exterior_forms, normal_vectors, jacobians,
+                           inverse_jacobians);
+  // quadrature_points, JxW_values,
+  // exterior_forms, normal_vectors);
+}
+
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::transform (
+  const VectorSlice<const std::vector<Tensor<1,dim> > > input,
+  VectorSlice<std::vector<Tensor<1,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+  AssertDimension (input.size(), output.size());
+
+  transform_fields(input, output, mapping_data, mapping_type);
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::transform (
+  const VectorSlice<const std::vector<DerivativeForm<1, dim ,spacedim>  > >  input,
+  VectorSlice<std::vector<Tensor<2,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+  AssertDimension (input.size(), output.size());
+
+  std::vector<DerivativeForm<1, spacedim,spacedim> > aux_output1(output.size());
+  VectorSlice< std::vector<DerivativeForm<1, spacedim,spacedim> > >  aux_output( aux_output1);
+
+  transform_differential_forms(input, aux_output, mapping_data, mapping_type);
+
+  for (unsigned int i=0; i<output.size(); i++)
+    output[i] = aux_output[i];
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void MappingFE<dim,spacedim,DH,VECTOR>::transform
+(const VectorSlice<const std::vector<Tensor<2, dim> > >     input,
+ VectorSlice<std::vector<Tensor<2,spacedim> > >             output,
+ const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+ const MappingType mapping_type) const
+{
+  AssertDimension (input.size(), output.size());
+
+  AssertThrow(false, ExcNotImplemented());
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+template < int rank >
+void MappingFE<dim,spacedim,DH,VECTOR>::transform_fields(
+  const VectorSlice<const std::vector<Tensor<rank,dim> > > input,
+  VectorSlice<std::vector<Tensor<rank,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_contravariant:
+    {
+      Assert (data.update_flags & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+      for (unsigned int i=0; i<output.size(); ++i)
+        output[i] = apply_transformation(data.contravariant[i], input[i]);
+
+      return;
+    }
+
+    case mapping_piola:
+    {
+      Assert (data.update_flags & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+      Assert (data.update_flags & update_volume_elements,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_volume_elements"));
+      Assert (rank==1, ExcMessage("Only for rank 1"));
+      for (unsigned int i=0; i<output.size(); ++i)
+        {
+          output[i] = apply_transformation(data.contravariant[i], input[i]);
+          output[i] /= data.volume_elements[i];
+        }
+      return;
+    }
+
+
+    //We still allow this operation as in the
+    //reference cell Derivatives are Tensor
+    //rather than DerivativeForm
+    case mapping_covariant:
+    {
+      Assert (data.update_flags & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+      for (unsigned int i=0; i<output.size(); ++i)
+        output[i] = apply_transformation(data.covariant[i], input[i]);
+
+      return;
+    }
+
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+template < int rank >
+void MappingFE<dim,spacedim,DH,VECTOR>::transform_differential_forms(
+  const VectorSlice<const std::vector<DerivativeForm<rank, dim,spacedim> > >    input,
+  VectorSlice<std::vector<DerivativeForm<rank, spacedim,spacedim> > > output,
+  const typename Mapping<dim,spacedim>::InternalDataBase &mapping_data,
+  const MappingType mapping_type) const
+{
+
+  AssertDimension (input.size(), output.size());
+  Assert (dynamic_cast<const InternalData *>(&mapping_data) != 0,
+          ExcInternalError());
+  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+
+  switch (mapping_type)
+    {
+    case mapping_covariant:
+    {
+      Assert (data.update_flags & update_contravariant_transformation,
+              typename FEValuesBase<dim>::ExcAccessToUninitializedField("update_contravariant_transformation"));
+
+      for (unsigned int i=0; i<output.size(); ++i)
+        output[i] = apply_transformation(data.covariant[i], input[i]);
+
+      return;
+    }
+    default:
+      Assert(false, ExcNotImplemented());
+    }
+
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+Point<spacedim>
+MappingFE<dim,spacedim,DH,VECTOR>::
+transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                             const Point<dim>                                 &p) const
+{
+//  Use the get_data function to create an InternalData with data vectors of
+//  the right size and transformation shape values already computed at point
+//  p. Make sure only one processor at a time performs this operations.
+
+  Threads::Mutex::ScopedLock lock(mutex);
+  update_internal_dofs(cell);
+  const Quadrature<dim> point_quadrature(p);
+  std::auto_ptr<InternalData>
+  mdata (dynamic_cast<InternalData *> (
+           get_data(update_transformation_values, point_quadrature)));
+
+  return this->transform_unit_to_real_cell_internal(*mdata);
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+Point<spacedim>
+MappingFE<dim,spacedim,DH,VECTOR>::
+transform_unit_to_real_cell_internal (const InternalData &data) const
+{
+  Point<spacedim> p_real;
+
+  for (unsigned int i=0; i<data.n_shape_functions; ++i)
+    {
+      unsigned int comp_i = fe->system_to_component_index(i).first;
+      if (fe_mask[comp_i])
+        p_real[fe_to_real[comp_i]] += local_dofs[i] * data.shape(0,i);
+    }
+
+  return p_real;
+}
+
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+Point<dim>
+MappingFE<dim,spacedim,DH,VECTOR>::
+transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+                             const Point<spacedim>                            &p) const
+{
+  {
+    Threads::Mutex::ScopedLock lock(mutex);
+    update_internal_dofs(cell);
+  }
+  // first a Newton iteration based on the real mapping. It uses the center
+  // point of the cell as a starting point
+  Point<dim> initial_p_unit;
+
+
+  for (unsigned int d=0; d<dim; ++d)
+    initial_p_unit[d] = 0.5;
+
+  // use the full mapping. in case the function above should have given us
+  // something back that lies outside the unit cell (that might happen
+  // because we may have given a point 'p' that lies inside the cell with
+  // the higher order mapping, but outside the Q1-mapped reference cell),
+  // then project it back into the reference cell in hopes that this gives
+  // a better starting point to the following iteration
+  initial_p_unit = GeometryInfo<dim>::project_to_unit_cell(initial_p_unit);
+
+  const Quadrature<dim> point_quadrature(initial_p_unit);
+
+  UpdateFlags update_flags = update_transformation_values|update_transformation_gradients;
+  if (spacedim>dim)
+    update_flags |= update_jacobian_grads;
+  std::auto_ptr<InternalData>
+  mdata (dynamic_cast<InternalData *> (
+           get_data(update_flags,point_quadrature)));
+
+  return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit, *mdata);
+
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+Point<dim>
+MappingFE<dim,spacedim,DH,VECTOR>::
+transform_real_to_unit_cell_internal
+(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+ const Point<spacedim>                            &p,
+ const Point<dim>                                 &initial_p_unit,
+ InternalData                                     &mdata) const
+{
+  {
+    Threads::Mutex::ScopedLock lock(mutex);
+
+    update_internal_dofs(cell);
+  }
+
+  const unsigned int n_shapes=mdata.shape_values.size();
+  Assert(n_shapes!=0, ExcInternalError());
+  AssertDimension (mdata.shape_derivatives.size(), n_shapes);
+
+
+  // Newton iteration to solve
+  // f(x)=p(x)-p=0
+  // x_{n+1}=x_n-[f'(x)]^{-1}f(x)
+
+  // The start value was set to be the
+  // linear approximation to the cell
+
+  // The shape values and derivatives
+  // of the mapping at this point are
+  // previously computed.
+
+  // For the <2,3>  case there is a
+  // template specialization.
+
+  // f(x)
+
+  //Point<spacedim> p_minus_F;
+  Point<dim> p_unit = initial_p_unit;
+  Point<dim> f;
+
+  compute_shapes_virtual(std::vector<Point<dim> > (1, p_unit), mdata);
+  Point<spacedim> p_real(transform_unit_to_real_cell_internal(mdata));
+
+  Tensor<1,spacedim> p_minus_F = p - p_real;
+
+  const double eps = 1.e-12*cell->diameter();
+  const unsigned int newton_iteration_limit = 20;
+
+  unsigned int newton_iteration=0;
+
+  while (p_minus_F.norm_square() > eps*eps)
+    {
+      // f'(x)
+      Point<spacedim>  DF[dim];
+      Tensor<2,dim>  df;
+
+      for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
+        {
+          const Tensor<1,dim> &grad_k = mdata.derivative(0,k);
+
+          unsigned int comp_k = fe->system_to_component_index(k).first;
+          if (fe_mask[comp_k])
+            for (unsigned int j=0; j<dim; ++j)
+              DF[j][fe_to_real[comp_k]] += local_dofs[k] * grad_k[j];
+        }
+
+      for (unsigned int j=0; j<dim; ++j)
+        {
+          f[j] = DF[j] * p_minus_F;
+          for (unsigned int l=0; l<dim; ++l)
+            df[j][l] = -DF[j] * DF[l];
+        }
+
+      // Solve  [f'(x)]d=f(x)
+      Tensor<1, dim> delta;
+      contract (delta, invert(df), static_cast<const Tensor<1,dim>&>(f));
+
+      // do a line search
+      double step_length = 1;
+      do
+        {
+          // update of p_unit. The
+          // spacedimth component of
+          // transformed point is simply
+          // ignored in codimension one
+          // case. When this component is
+          // not zero, then we are
+          // projecting the point to the
+          // surface or curve identified
+          // by the cell.
+          Point<dim> p_unit_trial = p_unit;
+          for (unsigned int i=0; i<dim; ++i)
+            p_unit_trial[i] -= step_length * delta[i];
+
+          // shape values and derivatives
+          // at new p_unit point
+          compute_shapes_virtual(std::vector<Point<dim> > (1, p_unit_trial), mdata);
+
+          // f(x)
+          Point<spacedim> p_real_trial = transform_unit_to_real_cell_internal(mdata);
+          //const Point<spacedim> f_trial = p - p_real_trial;
+          const Tensor<1,spacedim> f_trial = p - p_real_trial;
+
+          // see if we are making progress with the current step length
+          // and if not, reduce it by a factor of two and try again
+          if (f_trial.norm() < p_minus_F.norm())
+            {
+              p_real = p_real_trial;
+              p_unit = p_unit_trial;
+
+              p_minus_F = f_trial;
+              break;
+            }
+          else if (step_length > 0.05)
+            step_length /= 2;
+          else
+            {
+              std::cout << "Line search failed. With dim = " << dim << " spacedim = "
+                        << spacedim << std::endl;
+              goto failure;
+            }
+        }
+      while (true);
+
+      ++newton_iteration;
+      if (newton_iteration > newton_iteration_limit)
+        {
+          std::cout << "Too many newton iterations. With dim = " << dim << " spacedim = "
+                    << spacedim << std::endl;
+          goto failure;
+        }
+    }
+
+  return p_unit;
+
+  // if we get to the following label, then we have either run out
+  // of Newton iterations, or the line search has not converged.
+  // in either case, we need to give up, so throw an exception that
+  // can then be caught
+failure:
+  AssertThrow (false, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+
+  // ...the compiler wants us to return something, though we can
+  // of course never get here...
+  return Point<dim>();
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::compute_fill (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+  const unsigned int  n_q_points,
+  const typename QProjector<dim>::DataSetDescriptor  data_set,
+  const CellSimilarity::Similarity cell_similarity,
+  InternalData  &data,
+  std::vector<Point<spacedim> > &quadrature_points) const
+{
+  const UpdateFlags update_flags(data.current_update_flags());
+  Threads::Mutex::ScopedLock lock(mutex);
+  update_internal_dofs(cell);
+
+  // first compute quadrature points
+  if (update_flags & update_quadrature_points)
+    {
+      AssertDimension (quadrature_points.size(), n_q_points);
+
+      for (unsigned int point=0; point<n_q_points; ++point)
+        {
+          Point<spacedim> result;
+          const double *shape = &data.shape(point+data_set,0);
+
+          for (unsigned int k=0; k<data.n_shape_functions; ++k)
+            {
+              unsigned int comp_k = fe->system_to_component_index(k).first;
+              if (fe_mask[comp_k])
+                result[fe_to_real[comp_k]] += local_dofs[k] * shape[k];
+            }
+
+          quadrature_points[point] = result;
+        }
+    }
+
+
+  // then Jacobians
+  if (update_flags & update_contravariant_transformation)
+    {
+      AssertDimension (data.contravariant.size(), n_q_points);
+
+      // if the current cell is just a
+      // translation of the previous one, no
+      // need to recompute jacobians...
+      if (cell_similarity != CellSimilarity::translation)
+        {
+          std::fill(data.contravariant.begin(), data.contravariant.end(),
+                    DerivativeForm<1,dim,spacedim>());
+
+          Assert (data.n_shape_functions > 0, ExcInternalError());
+
+          for (unsigned int point=0; point<n_q_points; ++point)
+            {
+              const Tensor<1,dim> *data_derv =
+                &data.derivative(point+data_set, 0);
+
+              Tensor<1, dim> result[spacedim];
+
+              for (unsigned int k=0; k<data.n_shape_functions; ++k)
+                {
+                  unsigned int comp_k = fe->system_to_component_index(k).first;
+                  if (fe_mask[comp_k])
+                    result[fe_to_real[comp_k]] += local_dofs[k] * data_derv[k];
+                }
+
+              // write result into contravariant data. for
+              // j=dim in the case dim<spacedim, there will
+              // never be any nonzero data that arrives in
+              // here, so it is ok anyway because it was
+              // initialized to zero at the initialization
+              for (unsigned int i=0; i<spacedim; ++i)
+                {
+                  data.contravariant[point][i] = result[i];
+                }
+
+            }
+        }
+    }
+
+
+  if (update_flags & update_covariant_transformation)
+    {
+      AssertDimension (data.covariant.size(), n_q_points);
+      if (cell_similarity != CellSimilarity::translation)
+        for (unsigned int point=0; point<n_q_points; ++point)
+          data.covariant[point] = (data.contravariant[point]).covariant_form();
+    }
+
+  if (update_flags & update_volume_elements)
+    if (cell_similarity != CellSimilarity::translation)
+      for (unsigned int point=0; point<n_q_points; ++point)
+        data.volume_elements[point] = data.contravariant[point].determinant();
+
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::compute_fill_face (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell,
+  const unsigned int      face_no,
+  const unsigned int      subface_no,
+  const unsigned int      n_q_points,//npts
+  const typename QProjector<dim>::DataSetDescriptor data_set,
+  const std::vector<double>   &weights,
+  InternalData           &data,
+  std::vector<Point<spacedim> >    &quadrature_points,
+  std::vector<double>         &JxW_values,
+  std::vector<Tensor<1,spacedim> > &boundary_forms,
+  std::vector<Point<spacedim> > &normal_vectors,
+  std::vector<DerivativeForm<1,dim,spacedim> > &jacobians,
+  std::vector<DerivativeForm<1,spacedim,dim> > &inverse_jacobians) const
+{
+  compute_fill (cell, n_q_points, data_set, CellSimilarity::none,
+                data, quadrature_points);
+
+
+  const UpdateFlags update_flags(data.current_update_flags());
+
+  if (update_flags & update_boundary_forms)
+    {
+      AssertDimension (boundary_forms.size(), n_q_points);
+      if (update_flags & update_normal_vectors)
+        AssertDimension (normal_vectors.size(), n_q_points);
+      if (update_flags & update_JxW_values)
+        AssertDimension (JxW_values.size(), n_q_points);
+
+      // map the unit tangentials to the real cell. checking for d!=dim-1
+      // eliminates compiler warnings regarding unsigned int expressions <
+      // 0.
+      for (unsigned int d=0; d!=dim-1; ++d)
+        {
+          Assert (face_no+GeometryInfo<dim>::faces_per_cell*d <
+                  data.unit_tangentials.size(),
+                  ExcInternalError());
+          Assert (data.aux[d].size() <=
+                  data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d].size(),
+                  ExcInternalError());
+
+          transform (data.unit_tangentials[face_no+GeometryInfo<dim>::faces_per_cell*d],
+                     data.aux[d],
+                     data,
+                     mapping_contravariant);
+        }
+
+      // if dim==spacedim, we can use the unit tangentials to compute the
+      // boundary form by simply taking the cross product
+      if (dim == spacedim)
+        {
+          for (unsigned int i=0; i<n_q_points; ++i)
+            switch (dim)
+              {
+              case 1:
+                // in 1d, we don't have access to any of the data.aux
+                // fields (because it has only dim-1 components), but we
+                // can still compute the boundary form by simply
+                // looking at the number of the face
+                boundary_forms[i][0] = (face_no == 0 ?
+                                        -1 : +1);
+                break;
+              case 2:
+                cross_product (boundary_forms[i], data.aux[0][i]);
+                break;
+              case 3:
+                cross_product (boundary_forms[i], data.aux[0][i], data.aux[1][i]);
+                break;
+              default:
+                Assert(false, ExcNotImplemented());
+              }
+        }
+      else //(dim < spacedim)
+        {
+          // in the codim-one case, the boundary form results from the
+          // cross product of all the face tangential vectors and the cell
+          // normal vector
+          //
+          // to compute the cell normal, use the same method used in
+          // fill_fe_values for cells above
+          AssertDimension (data.contravariant.size(), n_q_points);
+
+          for (unsigned int point=0; point<n_q_points; ++point)
+            {
+              if (dim==1)
+                {
+                  // J is a tangent vector
+                  boundary_forms[point] = data.contravariant[point].transpose()[0];
+                  boundary_forms[point] /=
+                    (face_no == 0 ? -1. : +1.) * boundary_forms[point].norm();
+
+                }
+
+              if (dim==2)
+                {
+                  Tensor<1,spacedim> cell_normal;
+                  const DerivativeForm<1,spacedim,dim> DX_t =
+                    data.contravariant[point].transpose();
+                  cross_product(cell_normal,DX_t[0],DX_t[1]);
+                  cell_normal /= cell_normal.norm();
+
+                  // then compute the face normal from the face tangent
+                  // and the cell normal:
+                  cross_product (boundary_forms[point],
+                                 data.aux[0][point], cell_normal);
+
+                }
+
+            }
+        }
+
+
+
+      if (update_flags & (update_normal_vectors
+                          | update_JxW_values))
+        for (unsigned int i=0; i<boundary_forms.size(); ++i)
+          {
+            if (update_flags & update_JxW_values)
+              {
+                JxW_values[i] = boundary_forms[i].norm() * weights[i];
+
+                if (subface_no != numbers::invalid_unsigned_int)
+                  {
+                    const double area_ratio=GeometryInfo<dim>::subface_ratio(
+                                              cell->subface_case(face_no), subface_no);
+                    JxW_values[i] *= area_ratio;
+                  }
+              }
+
+            if (update_flags & update_normal_vectors)
+              normal_vectors[i] = Point<spacedim>(boundary_forms[i] / boundary_forms[i].norm());
+          }
+    }
+
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+unsigned int
+MappingFE<dim,spacedim,DH,VECTOR>::get_degree() const
+{
+  return fe->degree;
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+ComponentMask
+MappingFE<dim,spacedim,DH,VECTOR>::get_fe_mask() const
+{
+  return this->fe_mask;
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+Mapping<dim,spacedim> *
+MappingFE<dim,spacedim,DH,VECTOR>::clone () const
+{
+  return new MappingFE<dim,spacedim,DH,VECTOR>(*this);
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void
+MappingFE<dim,spacedim,DH,VECTOR>::update_internal_dofs (
+  const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+  if (euler_dof_handler == 0)
+    {
+      std::cout << "euler_dof_handler is empty!" << std::endl;
+      return;
+    }
+
+  typename DH::cell_iterator dof_cell(*cell, euler_dof_handler);
+  Assert (dof_cell->active() == true, ExcInactiveCell());
+
+  dof_cell->get_dof_indices(dof_indices);
+
+  for (unsigned int i=0; i<local_dofs.size(); ++i)
+    {
+      local_dofs[i] = (*euler_vector)(dof_indices[i]);
+    }
+}
+
+
+template<int dim, int spacedim, class DH, class VECTOR>
+void MappingFE<dim,spacedim,DH,VECTOR>::update_euler_vector_using_triangulation
+(VECTOR &vector)
+{
+  if ( fe->has_support_points() )
+    {
+      std::vector<Point<dim> > support_points = fe->get_unit_support_points();
+      typename DH::active_cell_iterator cell;
+      Quadrature<dim> quad(support_points);
+
+      MappingQ<dim,spacedim> map_q(fe->degree);
+      FEValues<dim,spacedim> fe_v(map_q, *fe, quad, update_quadrature_points);
+      std::vector<unsigned int> dofs(fe->dofs_per_cell);
+
+      AssertDimension(fe->dofs_per_cell, support_points.size());
+      Assert(fe->is_primitive(), ExcMessage("FE is not Primitive! This won't work."));
+
+      for (cell = euler_dof_handler->begin_active(); cell != euler_dof_handler->end(); ++cell)
+        {
+          fe_v.reinit(cell);
+          cell->get_dof_indices(dofs);
+          const std::vector<Point<spacedim> > &points = fe_v.get_quadrature_points();
+          for (unsigned int q = 0; q < points.size(); ++q)
+            {
+              unsigned int comp = fe->system_to_component_index(q).first;
+              vector(dofs[q]) = points[q][comp];
+            }
+        }
+
+    }
+  else
+    {
+      // Construct a MappingFE with an FEQ
+      FESystem<dim,spacedim> feq(FE_Q<dim,spacedim>(fe->degree), spacedim);
+      DH dhq(euler_dof_handler->get_tria());
+      dhq.distribute_dofs(feq);
+      VECTOR eulerq(dhq.n_dofs());
+      const ComponentMask maskq(spacedim, true);
+      MappingFE<dim,spacedim,DH,VECTOR> newfe(eulerq, dhq, maskq);
+
+      newfe.update_euler_vector_using_triangulation(eulerq);
+
+      FullMatrix<double> transfer(fe->dofs_per_cell, feq.dofs_per_cell);
+      std::vector<Point<dim> > points = feq.get_unit_support_points();
+
+      // Here construct the matrix!!!!
+      for (unsigned int i=0; i<fe->dofs_per_cell; ++i)
+        {
+          for (unsigned int j=0; j<points.size(); ++j)
+            {
+              if (fe->system_to_component_index(i).first
+                  ==
+                  feq.system_to_component_index(j).first)
+                transfer(j,i) = fe->shape_value(i, points[j]);
+            }
+        }
+      VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector);
+    }
+}
+
+
+
+
+// explicit instantiations
+#include "mapping_fe.inst"
+
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/fe/mapping_fe.inst.in b/source/fe/mapping_fe.inst.in
new file mode 100644 (file)
index 0000000..a06af31
--- /dev/null
@@ -0,0 +1,26 @@
+// ---------------------------------------------------------------------
+// $Id: mapping_fe.inst.in 30049 2013-07-18 19:42:40Z maier $
+//
+// Copyright (C) 1998 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+  {
+#if deal_II_dimension <= deal_II_space_dimension       
+    template class MappingFE<deal_II_dimension, deal_II_space_dimension, dealii::DoFHandler<deal_II_dimension, deal_II_space_dimension>, dealii::Vector<double> >;
+#endif
+  }
+
diff --git a/tests/fe/mapping_fe_real_to_unit_q1.cc b/tests/fe/mapping_fe_real_to_unit_q1.cc
new file mode 100644 (file)
index 0000000..10097bb
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// on a somewhat deformed cube, verify that if we push forward a bunch
+// of points from the reference to the real cell and then call
+// MappingFE::transform_unit_to_real_cell that we get the same point as
+// we had in the beginning.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/component_mask.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/fe_system.h>
+
+using namespace dealii;
+
+
+template<int dim, int spacedim>
+void test_real_to_unit_cell()
+{
+       deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl;
+
+ Triangulation<dim, spacedim>   triangulation;
+ GridGenerator::hyper_cube (triangulation);
+
+ const unsigned int n_points = 4;
+ std::vector< Point<dim> > unit_points(Utilities::fixed_power<dim>(n_points));
+
+ switch (dim)
+   {
+   case 1:
+     for (unsigned int x=0; x<n_points; ++x)
+       unit_points[x][0] = double(x)/double(n_points);
+     break;
+
+   case 2:
+     for (unsigned int x=0; x<n_points; ++x)
+       for (unsigned int y=0; y<n_points; ++y)
+         {
+           unit_points[y * n_points + x][0] = double(x)/double(n_points);
+           unit_points[y * n_points + x][1] = double(y)/double(n_points);
+         }
+     break;
+
+   case 3:
+     for (unsigned int x=0; x<n_points; ++x)
+       for (unsigned int y=0; y<n_points; ++y)
+         for (unsigned int z=0; z<n_points; ++z)
+           {
+             unit_points[z * n_points * n_points + y * n_points + x][0] = double(x)/double(n_points);
+             unit_points[z * n_points * n_points + y * n_points + x][1] = double(y)/double(n_points);
+             unit_points[z * n_points * n_points + y * n_points + x][2] = double(z)/double(n_points);
+           }
+     break;
+   }
+
+ const FE_Q<dim,spacedim> feq(1);
+ const FESystem<dim,spacedim> fesystem(feq, spacedim);
+ DoFHandler<dim,spacedim> dhq(triangulation);
+ dhq.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhq.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+
+ map.update_euler_vector_using_triangulation(eulerq);
+
+ typename Triangulation<dim, spacedim >::active_cell_iterator
+ cell = triangulation.begin_active();
+
+  //Move a vertex a little bit
+  const unsigned int n_dx = 5;
+  const double dx = 0.4/n_dx;
+  Point<spacedim> direction;
+  for (unsigned int j=0; j<spacedim; ++j)
+    direction[j]=dx;
+
+  // in a loop, make the cell more
+  // and more distorted
+  for (unsigned int j=0; j<n_dx; ++j)
+  {
+    deallog << "Vertex displacement: " << double(j)*direction <<  std::endl;
+    cell->vertex(0) = double(j)*direction;
+
+    for (unsigned int i=0; i<unit_points.size(); ++i)
+    {
+         // for each of the points,
+         // verify that if we apply
+         // the forward map and then
+         // pull back that we get
+         // the same point again
+     
+      const Point<spacedim> p = map.transform_unit_to_real_cell(cell,unit_points[i]);
+      const Point<dim> p_unit = map.transform_real_to_unit_cell(cell,p);
+     
+      Assert (unit_points[i].distance(p_unit) < 1e-10, ExcInternalError());
+    }
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+       std::ofstream logfile ("output");
+       deallog.attach(logfile);
+       deallog.depth_console(0);
+       deallog.threshold_double(1.e-10);
+
+       test_real_to_unit_cell<1,1>();
+       test_real_to_unit_cell<2,2>();
+       test_real_to_unit_cell<3,3>();
+
+       test_real_to_unit_cell<1,2>();
+       test_real_to_unit_cell<1,3>();
+       test_real_to_unit_cell<2,3>();
+
+       return 0;
+}
diff --git a/tests/fe/mapping_fe_real_to_unit_q1.output b/tests/fe/mapping_fe_real_to_unit_q1.output
new file mode 100644 (file)
index 0000000..c68c1fb
--- /dev/null
@@ -0,0 +1,43 @@
+
+DEAL::dim=1, spacedim=1
+DEAL::Vertex displacement: 0.00000
+DEAL::Vertex displacement: 0.0800000
+DEAL::Vertex displacement: 0.160000
+DEAL::Vertex displacement: 0.240000
+DEAL::Vertex displacement: 0.320000
+DEAL::OK
+DEAL::dim=2, spacedim=2
+DEAL::Vertex displacement: 0.00000 0.00000
+DEAL::Vertex displacement: 0.0800000 0.0800000
+DEAL::Vertex displacement: 0.160000 0.160000
+DEAL::Vertex displacement: 0.240000 0.240000
+DEAL::Vertex displacement: 0.320000 0.320000
+DEAL::OK
+DEAL::dim=3, spacedim=3
+DEAL::Vertex displacement: 0.00000 0.00000 0.00000
+DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000
+DEAL::Vertex displacement: 0.160000 0.160000 0.160000
+DEAL::Vertex displacement: 0.240000 0.240000 0.240000
+DEAL::Vertex displacement: 0.320000 0.320000 0.320000
+DEAL::OK
+DEAL::dim=1, spacedim=2
+DEAL::Vertex displacement: 0.00000 0.00000
+DEAL::Vertex displacement: 0.0800000 0.0800000
+DEAL::Vertex displacement: 0.160000 0.160000
+DEAL::Vertex displacement: 0.240000 0.240000
+DEAL::Vertex displacement: 0.320000 0.320000
+DEAL::OK
+DEAL::dim=1, spacedim=3
+DEAL::Vertex displacement: 0.00000 0.00000 0.00000
+DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000
+DEAL::Vertex displacement: 0.160000 0.160000 0.160000
+DEAL::Vertex displacement: 0.240000 0.240000 0.240000
+DEAL::Vertex displacement: 0.320000 0.320000 0.320000
+DEAL::OK
+DEAL::dim=2, spacedim=3
+DEAL::Vertex displacement: 0.00000 0.00000 0.00000
+DEAL::Vertex displacement: 0.0800000 0.0800000 0.0800000
+DEAL::Vertex displacement: 0.160000 0.160000 0.160000
+DEAL::Vertex displacement: 0.240000 0.240000 0.240000
+DEAL::Vertex displacement: 0.320000 0.320000 0.320000
+DEAL::OK
diff --git a/tests/fe/mapping_fe_real_to_unit_q5_curved.cc b/tests/fe/mapping_fe_real_to_unit_q5_curved.cc
new file mode 100644 (file)
index 0000000..486c306
--- /dev/null
@@ -0,0 +1,151 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2006 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// on a somewhat deformed cube, verify that if we push forward a bunch
+// of points from the reference to the real cell and then call
+// MappingFE::transform_unit_to_real_cell that we get the same point as
+// we had in the beginning.
+
+// We use a Q5 mapping but this time we
+// actually curve one boundary of the cell which ensures that the
+// mapping is really higher order than just Q1
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/component_mask.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <fstream>
+
+using namespace dealii;
+
+
+template<int dim, int spacedim>
+void test_real_to_unit_cell()
+{
+       deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl;
+
+  // define a boundary that fits the
+  // the vertices of the hyper cube
+  // we're going to create below
+  HyperBallBoundary<dim,spacedim> boundary (Point<spacedim>(),
+                                            std::sqrt(1.*dim));
+
+  Triangulation<dim, spacedim>   triangulation;
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+
+  // set the boundary indicator for
+  // one face of the single cell
+  triangulation.set_boundary (1, boundary);
+  triangulation.begin_active()->face(0)->set_boundary_indicator (1);
+
+
+ const unsigned int n_points = 5;
+ std::vector< Point<dim> > unit_points(Utilities::fixed_power<dim>(n_points));
+
+ switch (dim)
+   {
+   case 1:
+     for (unsigned int x=0; x<n_points; ++x)
+       unit_points[x][0] = double(x)/double(n_points);
+     break;
+
+   case 2:
+     for (unsigned int x=0; x<n_points; ++x)
+       for (unsigned int y=0; y<n_points; ++y)
+         {
+           unit_points[y * n_points + x][0] = double(x)/double(n_points);
+           unit_points[y * n_points + x][1] = double(y)/double(n_points);
+         }
+     break;
+
+   case 3:
+     for (unsigned int x=0; x<n_points; ++x)
+       for (unsigned int y=0; y<n_points; ++y)
+         for (unsigned int z=0; z<n_points; ++z)
+           {
+             unit_points[z * n_points * n_points + y * n_points + x][0] = double(x)/double(n_points);
+             unit_points[z * n_points * n_points + y * n_points + x][1] = double(y)/double(n_points);
+             unit_points[z * n_points * n_points + y * n_points + x][2] = double(z)/double(n_points);
+           }
+     break;
+   }
+
+ const FE_Q<dim,spacedim> feq(5);
+ const FESystem<dim,spacedim> fesystem(feq, spacedim);
+ DoFHandler<dim,spacedim> dhq(triangulation);
+ dhq.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhq.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+
+ map.update_euler_vector_using_triangulation(eulerq);
+
+ typename Triangulation<dim, spacedim >::active_cell_iterator
+ cell = triangulation.begin_active();
+
+      for (unsigned int i=0; i<unit_points.size(); ++i)
+      {
+                 // for each of the points,
+                 // verify that if we apply
+                 // the forward map and then
+                 // pull back that we get
+                 // the same point again
+
+       const Point<spacedim> p = map.transform_unit_to_real_cell(cell,unit_points[i]);
+
+       const Point<dim> p_unit = map.transform_real_to_unit_cell(cell,p);
+
+       AssertThrow (unit_points[i].distance(p_unit) < 1e-10, ExcInternalError());
+     }
+
+  deallog << "OK" << std::endl;
+
+}
+
+
+int
+main()
+{
+       std::ofstream logfile ("output");
+       deallog.attach(logfile);
+       deallog.depth_console(0);
+       deallog.threshold_double(1.e-10);
+
+       test_real_to_unit_cell<1,1>();
+       test_real_to_unit_cell<2,2>();
+  test_real_to_unit_cell<3,3>();
+
+       test_real_to_unit_cell<1,2>();
+       test_real_to_unit_cell<2,3>();
+
+
+  //test_real_to_unit_cell<1,3>();
+       return 0;
+}
diff --git a/tests/fe/mapping_fe_real_to_unit_q5_curved.output b/tests/fe/mapping_fe_real_to_unit_q5_curved.output
new file mode 100644 (file)
index 0000000..5693e8c
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::dim=1, spacedim=1
+DEAL::OK
+DEAL::dim=2, spacedim=2
+DEAL::OK
+DEAL::dim=3, spacedim=3
+DEAL::OK
+DEAL::dim=1, spacedim=2
+DEAL::OK
+DEAL::dim=2, spacedim=3
+DEAL::OK

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.