]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Renamed MappingFE to MappingFEField
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 6 Apr 2015 18:06:32 +0000 (20:06 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 16 Apr 2015 14:20:12 +0000 (16:20 +0200)
Fixed documentation.

Fixed indentation.

Inverted DH and VECTOR in template declaration.

Renamed get_fe_mask to get_component_mask

Added entry in changes.h

fixed: destructor, useless lines, commented code. Renamed tests

added missing documentation and solved the thread safe issue

eliminated if statement in compute_shapes_virtual

solved all the problems highlighted in https://github.com/dealii/dealii/pull/727

const quadrature points in update_euler_vector_using_triangulation

Added two test with a non trivial component mask (mapping_fe_field_real_to_unit_q1_mask.cc) and fixed if statement in update_euler_vector_using_triangulation

Added assert dimension in update_euler_vector_using_triangulation

start working on else statement in update_euler

added one test with fe_bernstein (still failing)

Added tests from q1 to q5 and from b1 to b5

Fixed doxygen comments.

Added a few comments and fixed a const reference.

fixed inital guess in transform_real_to_unit

added hyperballboundary in test b2_mask

29 files changed:
doc/news/changes.h
include/deal.II/fe/mapping_fe_field.h [moved from include/deal.II/fe/mapping_fe.h with 72% similarity]
source/fe/CMakeLists.txt
source/fe/mapping_fe_field.cc [moved from source/fe/mapping_fe.cc with 80% similarity]
source/fe/mapping_fe_field.inst.in [moved from source/fe/mapping_fe.inst.in with 82% similarity]
tests/fe/mapping_fe_field_real_to_unit_b1.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b1.output [moved from tests/fe/mapping_fe_real_to_unit_q1.output with 100% similarity]
tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b2_curved.output [moved from tests/fe/mapping_fe_real_to_unit_q5_curved.output with 100% similarity]
tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b2_mask.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b3_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b4_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_b5_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q1.cc [moved from tests/fe/mapping_fe_real_to_unit_q1.cc with 95% similarity]
tests/fe/mapping_fe_field_real_to_unit_q1.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q2_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q2_mask.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q3_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q4_curved.output [new file with mode: 0644]
tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc [moved from tests/fe/mapping_fe_real_to_unit_q5_curved.cc with 95% similarity]
tests/fe/mapping_fe_field_real_to_unit_q5_curved.output [new file with mode: 0644]

index aae3579f89797b8a4341543ac3978a40324aa0ba..2c3d4c119bf8efd3efd0b279063029ddaf0d7be0 100644 (file)
@@ -511,6 +511,12 @@ inconvenience this causes.
   (Wolfgang Bangerth, 2015/04/08)
   </li>
 
+  <li> New: A new MappingFEField() class has been added to the library
+  that generalizes MappingQEulerian to allow arbitrary FiniteElements.
+  <br>
+  (Marco Tezzele, Luca Heltai, 2015/04/06)
+  <li>
+
   <li> Changed: In the spirit of the changes made to the distinction
   between Point and Tensor objects discussed above, the first argument
   to GridTools::shift() has been changed from a Point to a Tensor@<1,dim@>.
similarity index 72%
rename from include/deal.II/fe/mapping_fe.h
rename to include/deal.II/fe/mapping_fe_field.h
index 308c0a82dbe4b6ca47d1e19c80bdff32df27e663..7e065af27ca906334464a8e26517b50054cf10f9 100644 (file)
@@ -32,29 +32,29 @@ DEAL_II_NAMESPACE_OPEN
 /*@{*/
 
 /**
- * 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 displacements.
- * 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).
+ * The MappingFEField is a generalization of the MappingQEulerian
+ * class, for arbitrary vector finite elements. The two main
+ * differences are that this class uses a vector of absolute positions
+ * instead of a vector of displacements, and it allows for arbitrary
+ * FiniteElement types, instead of only FE_Q.
  *
+ * This class effectively decouples the topology from the geometry, by
+ * relegating all geometrical information to some components of a
+ * FiniteElement vector field. The components that are used for the
+ * geometry can be arbitrarily selected at construction time.
+ *
+ * The idea is to consider the Triangulation as a parameter
+ * configuration space, on which we  construct an arbitrary
+ * geometrical mapping, using the instruments of the deal.II library:
+ * a vector of degrees of freedom, a DoFHandler associated to the
+ * geometry of the problem and a ComponentMask that tells us which
+ * components of the FiniteElement to use for the mapping.
+ *
+ * Typically, the DoFHandler operates on a finite element that is
+ * constructed as a system element (FESystem) from continuous FE_Q()
+ * (for iso-parametric discretizations) or FE_Bernstein() (for
+ * iso-geometric discretizations) objects. An example is shown below:
  *
- * 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);
@@ -62,50 +62,86 @@ DEAL_II_NAMESPACE_OPEN
  *    dhq.distribute_dofs(fesystem);
  *    Vector<double> eulerq(dhq.n_dofs());
  *    const ComponentMask mask(spacedim, true);
- *    MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+ *    MappingFEField<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>
+          class VECTOR=Vector<double>,
+          class DH=DoFHandler<dim,spacedim> >
+class MappingFEField : 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.
+   * configuration.
+   *
+   * In general this class decouples geometry from topology, allowing
+   * users to define geometries which are only topologically
+   * equivalent to the underlying Triangulation, but which may
+   * otherwise be arbitrary. Differently from what happens in
+   * MappingQEulerian, the FiniteElement field which is passed to the
+   * constructor is interpreted as an absolute geometrical
+   * configuration, therefore one has to make sure that the
+   * euler_vector actually represents a valid geometry (i.e., one with
+   * no inverted cells, or with no zero-volume cells).
+   *
+   * In order to facilitate the user, we provide an
+   * update_euler_vector_using_triangulation() method to initialize a
+   * euler_vector in a way similar to what happens in
+   * MappingQEulerian, by creating a euler_vector which interpolates
+   * the underlying Triangulation.
+   *
+   * If the underlying FiniteElement is a system of FE_Q(), and
+   * euler_vector is initialized using
+   * update_euler_vector_using_triangulation(), then this class is in
+   * all respects identical to MappingQ().
+   *
+   * The optional ComponentMask argument can be used to specify what
+   * components of the FiniteElement to use for the geometrical
+   * transformation. If no mask is specified at construction time,
+   * then a default one is used, which makes this class works in the
+   * same way of MappingQEulerian(), i.e., the first spacedim
+   * components of the FiniteElement are assumed to represent the
+   * geometry of the problem.
+   *
+   * Notice that if a mask is specified, it has to match in size the
+   * underlying FiniteElement, and it has to have exactly spacedim
+   * non-zero elements, indicating the components (in order) of the
+   * FiniteElement which will be used for the geometry.
+   *
+   * If an incompatible mask is passed, an exception is thrown.
    */
-  MappingFE (const VECTOR  &euler_vector,
-             const DH      &euler_dof_handler,
-             const ComponentMask mask=ComponentMask());
+  MappingFEField (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);
+  MappingFEField (const MappingFEField<dim,spacedim,VECTOR,DH> &mapping);
 
   /**
-   * Destructor.
+   * Helper function to fill the given euler vector with information
+   * coming from the triangulation, by interpolating the geometry of
+   * the Triangulation associated with the DoFHandler used at
+   * construction time.
+   *
+   * The resulting map is guaranteed to be interpolatory at the
+   * vertices of the Triangulation. Notice that this may or may not be
+   * meaningful, depending on the FiniteElement you have used to
+   * construct this MappingFEField.
+   *
+   * If the underlying FiniteElement is a system of FE_Q(), and
+   * euler_vector is initialized using this function, then this class
+   * is in all respects identical to MappingQ().
    */
-  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);
+  void update_euler_vector_using_triangulation(VECTOR &vector) const;
 
 
 
@@ -132,19 +168,27 @@ public:
     const typename Triangulation<dim,spacedim>::cell_iterator &cell,
     const Point<spacedim>                            &p) const;
 
-
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   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;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   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;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   virtual
   void
   transform (const VectorSlice<const std::vector<Tensor<2, dim> > >     input,
@@ -164,7 +208,7 @@ public:
    * Return the ComponentMask of the mapping, i.e. which components to use for
    * the mapping.
    */
-  ComponentMask get_fe_mask () const;
+  ComponentMask get_component_mask () const;
 
   /**
    * Return a pointer to a copy of the present object. The caller of this copy
@@ -336,20 +380,26 @@ public:
     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.
+     * 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;
 
+    /**
+     * Stores the mask given at construction time. If no mask was
+     * specified at construction time, then a default one is used,
+     * which makes this class works in the same way of
+     * MappingQEulerian(), i.e., the first spacedim components of the
+     * FiniteElement are used for the euler_vector and the euler_dh.
+     *
+     * If a mask is specified, then it has to match the underlying
+     * FiniteElement, and it has to have exactly spacedim non-zero
+     * elements, indicating the components (in order) of the
+     * FiniteElement which will be used for the euler vector and the
+     * euler dof handler.
+     */
     ComponentMask mask;
   };
 
@@ -537,7 +587,7 @@ protected:
   * Reference to the vector of shifts.
   */
 
-  SmartPointer<const VECTOR, MappingFE<dim,spacedim,DH,VECTOR> >euler_vector;
+  SmartPointer<const VECTOR, MappingFEField<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
@@ -545,13 +595,13 @@ protected:
    * 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;
+  SmartPointer<const FiniteElement<dim,spacedim>, MappingFEField<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;
+  SmartPointer<const DH,MappingFEField<dim,spacedim,DH,VECTOR> >euler_dof_handler;
 
 
 
@@ -563,49 +613,78 @@ private:
   void update_internal_dofs(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
 
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   mutable std::vector<double> local_dofs;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   mutable std::vector<types::global_dof_index> dof_indices;
 
   /**
    * Mutex to protect local_dofs.
    */
-
   mutable Threads::Mutex mutex;
 
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   virtual void
   compute_shapes_virtual (const std::vector<Point<dim> > &unit_points,
-                          typename MappingFE<dim, spacedim>::InternalData &data) const;
+                          typename MappingFEField<dim, spacedim>::InternalData &data) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   UpdateFlags
   update_once (const UpdateFlags in) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   UpdateFlags
   update_each (const UpdateFlags in) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   void
   compute_data (const UpdateFlags      update_flags,
                 const Quadrature<dim>  &q,
                 const unsigned int     n_original_q_points,
                 InternalData           &data) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   void
   compute_face_data (const UpdateFlags      update_flags,
                      const Quadrature<dim>  &q,
                      const unsigned int     n_original_q_points,
                      InternalData           &data) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   virtual
   typename Mapping<dim,spacedim>::InternalDataBase *
   get_data (const UpdateFlags,
             const Quadrature<dim> &quadrature) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   virtual
   typename Mapping<dim,spacedim>::InternalDataBase *
   get_face_data (const UpdateFlags flags,
                  const Quadrature<dim-1>& quadrature) const;
 
+  /** Reimplemented from Mapping. See the documentation of the base class for
+   * detailed information.
+   */
   virtual
   typename Mapping<dim,spacedim>::InternalDataBase *
   get_subface_data (const UpdateFlags flags,
@@ -630,9 +709,9 @@ private:
 
 
   /**
-   * Declare other MappingFE classes friends.
+   * Declare other MappingFEField classes friends.
    */
-  template <int,int,class,class> friend class MappingFE;
+  template <int,int,class,class> friend class MappingFEField;
 };
 
 /*@}*/
@@ -645,8 +724,8 @@ private:
 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
+MappingFEField<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,
@@ -659,8 +738,8 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::shape (const unsigned int qpoin
 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)
+MappingFEField<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,
@@ -672,7 +751,7 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::shape (const unsigned int qpoin
 template<int dim, int spacedim, class DH, class VECTOR>
 inline
 Tensor<1,dim>
-MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int qpoint,
+MappingFEField<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(),
@@ -686,7 +765,7 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int
 template<int dim, int spacedim, class DH, class VECTOR>
 inline
 Tensor<1,dim> &
-MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int qpoint,
+MappingFEField<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(),
@@ -699,7 +778,7 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::derivative (const unsigned int
 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,
+MappingFEField<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(),
@@ -713,7 +792,7 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::second_derivative (const unsign
 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,
+MappingFEField<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(),
@@ -726,7 +805,7 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::second_derivative (const unsign
 template <int dim, int spacedim, class DH, class VECTOR>
 inline
 bool
-MappingFE<dim,spacedim,DH,VECTOR>::preserves_vertex_locations () const
+MappingFEField<dim,spacedim,DH,VECTOR>::preserves_vertex_locations () const
 {
   return false;
 }
index def0a6685aa555b552a8d726550ae89257f3c3e2..9d7c756970b478f10134988208fc741de8ea8120 100644 (file)
@@ -49,7 +49,7 @@ SET(_src
   mapping_c1.cc
   mapping_cartesian.cc
   mapping.cc
-  mapping_fe.cc
+  mapping_fe_field.cc
   mapping_q1.cc
   mapping_q1_eulerian.cc
   mapping_q.cc
@@ -90,7 +90,7 @@ SET(_inst
   mapping_c1.inst.in
   mapping_cartesian.inst.in
   mapping.inst.in
-  mapping_fe.inst.in
+  mapping_fe_field.inst.in
   mapping_q1_eulerian.inst.in
   mapping_q1.inst.in
   mapping_q_eulerian.inst.in
similarity index 80%
rename from source/fe/mapping_fe.cc
rename to source/fe/mapping_fe_field.cc
index c8343e5f0323c5b6ed8d4b7a1d1d8759c81b2ce1..c2a199fc73c24c2f7f4b4ef9f8eeb0181b414648 100644 (file)
@@ -27,7 +27,7 @@
 #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/mapping_fe_field.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/mapping.h>
 #include <deal.II/fe/mapping_q1.h>
@@ -44,8 +44,8 @@
 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,
+template<int dim, int spacedim, class VECTOR, class DH>
+MappingFEField<dim,spacedim,VECTOR,DH>::InternalData::InternalData (const FiniteElement<dim,spacedim> &fe,
     const ComponentMask mask)
   :
   n_shape_functions (fe.dofs_per_cell),
@@ -54,18 +54,18 @@ MappingFE<dim,spacedim,DH,VECTOR>::InternalData::InternalData (const FiniteEleme
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 std::size_t
-MappingFE<dim,spacedim,DH,VECTOR>::InternalData::memory_consumption () const
+MappingFEField<dim,spacedim,VECTOR,DH>::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)
+template<int dim, int spacedim, class VECTOR, class DH>
+MappingFEField<dim,spacedim,VECTOR,DH>::MappingFEField (const VECTOR  &euler_vector,
+                                                        const DH      &euler_dof_handler,
+                                                        const ComponentMask mask)
   :
   euler_vector(&euler_vector),
   fe(&euler_dof_handler.get_fe()),
@@ -76,7 +76,6 @@ MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const VECTOR  &euler_vector,
           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)
     {
@@ -87,8 +86,8 @@ MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const VECTOR  &euler_vector,
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
-MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const MappingFE<dim,spacedim,DH,VECTOR> &mapping)
+template<int dim, int spacedim, class VECTOR, class DH>
+MappingFEField<dim,spacedim,VECTOR,DH>::MappingFEField (const MappingFEField<dim,spacedim,VECTOR,DH> &mapping)
   :
   euler_vector(mapping.euler_vector),
   fe(mapping.fe),
@@ -100,44 +99,34 @@ MappingFE<dim,spacedim,DH,VECTOR>::MappingFE (const MappingFE<dim,spacedim,DH,VE
 {}
 
 
-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>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::compute_shapes_virtual (
+MappingFEField<dim,spacedim,VECTOR,DH>::compute_shapes_virtual (
   const std::vector<Point<dim> > &unit_points,
-  typename MappingFE<dim, spacedim>::InternalData &data) const
+  typename MappingFEField<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]);
-      }
+  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>
+template<int dim, int spacedim, class VECTOR, class DH>
 UpdateFlags
-MappingFE<dim,spacedim,DH,VECTOR>::update_once (const UpdateFlags in) const
+MappingFEField<dim,spacedim,VECTOR,DH>::update_once (const UpdateFlags in) const
 {
   UpdateFlags out = UpdateFlags(in & (update_transformation_values
                                       | update_transformation_gradients));
@@ -162,9 +151,9 @@ MappingFE<dim,spacedim,DH,VECTOR>::update_once (const UpdateFlags in) const
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 UpdateFlags
-MappingFE<dim,spacedim,DH,VECTOR>::update_each (const UpdateFlags in) const
+MappingFEField<dim,spacedim,VECTOR,DH>::update_each (const UpdateFlags in) const
 {
   // Select flags of concern for the
   // transformation.
@@ -187,11 +176,11 @@ MappingFE<dim,spacedim,DH,VECTOR>::update_each (const UpdateFlags in) const
   // set operations. this leads to
   // some circular logic. the only
   // way to treat this is to
-  // iterate. since there are 4
+  // iterate. since there are 5
   // if-clauses in the loop, it will
-  // take at most 3 iterations to
+  // take at most 4 iterations to
   // converge. do them:
-  for (unsigned int i=0; i<4; ++i)
+  for (unsigned int i=0; i<5; ++i)
     {
       // The following is a little incorrect:
       // If not applied on a face,
@@ -236,12 +225,12 @@ MappingFE<dim,spacedim,DH,VECTOR>::update_each (const UpdateFlags in) const
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 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
+MappingFEField<dim,spacedim,VECTOR,DH>::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();
 
@@ -270,16 +259,15 @@ MappingFE<dim,spacedim,DH,VECTOR>::compute_data (const UpdateFlags      update_f
     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>
+template<int dim, int spacedim, class VECTOR, class DH>
 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
+MappingFEField<dim,spacedim,VECTOR,DH>::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);
 
@@ -341,10 +329,10 @@ MappingFE<dim,spacedim,DH,VECTOR>::compute_face_data (const UpdateFlags update_f
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 typename Mapping<dim,spacedim>::InternalDataBase *
-MappingFE<dim,spacedim,DH,VECTOR>::get_data (const UpdateFlags update_flags,
-                                             const Quadrature<dim> &quadrature) const
+MappingFEField<dim,spacedim,VECTOR,DH>::get_data (const UpdateFlags update_flags,
+                                                  const Quadrature<dim> &quadrature) const
 {
   InternalData *data = new InternalData(*fe, fe_mask);
   this->compute_data (update_flags, quadrature,
@@ -354,10 +342,10 @@ MappingFE<dim,spacedim,DH,VECTOR>::get_data (const UpdateFlags update_flags,
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 typename Mapping<dim,spacedim>::InternalDataBase *
-MappingFE<dim,spacedim,DH,VECTOR>::get_face_data (const UpdateFlags update_flags,
-                                                  const Quadrature<dim-1>& quadrature) const
+MappingFEField<dim,spacedim,VECTOR,DH>::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));
@@ -368,10 +356,10 @@ MappingFE<dim,spacedim,DH,VECTOR>::get_face_data (const UpdateFlags update_flags
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 typename Mapping<dim,spacedim>::InternalDataBase *
-MappingFE<dim,spacedim,DH,VECTOR>::get_subface_data (const UpdateFlags update_flags,
-                                                     const Quadrature<dim-1>& quadrature) const
+MappingFEField<dim,spacedim,VECTOR,DH>::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));
@@ -382,11 +370,11 @@ MappingFE<dim,spacedim,DH,VECTOR>::get_subface_data (const UpdateFlags update_fl
 }
 
 
-// Note that the CellSimilarity flag is modifyable, since MappingFE can need to
+// Note that the CellSimilarity flag is modifyable, since MappingFEField can need to
 // recalculate data even when cells are similar.
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_values (
+MappingFEField<dim,spacedim,VECTOR,DH>::fill_fe_values (
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const Quadrature<dim>                                     &q,
   typename Mapping<dim,spacedim>::InternalDataBase          &mapping_data,
@@ -406,14 +394,6 @@ MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_values (
   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;
 
@@ -567,9 +547,9 @@ MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_values (
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_face_values (
+MappingFEField<dim,spacedim,VECTOR,DH>::fill_fe_face_values (
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const unsigned int       face_no,
   const Quadrature<dim-1> &q,
@@ -580,16 +560,9 @@ MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_face_values (
   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);
@@ -608,30 +581,23 @@ MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_face_values (
                            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>
+template<int dim, int spacedim, class VECTOR, class DH>
 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
+MappingFEField<dim,spacedim,VECTOR,DH>::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
 {
-  //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,
@@ -653,16 +619,13 @@ MappingFE<dim,spacedim,DH,VECTOR>::fill_fe_subface_values (const typename Triang
                            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>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::transform (
+MappingFEField<dim,spacedim,VECTOR,DH>::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,
@@ -675,9 +638,9 @@ MappingFE<dim,spacedim,DH,VECTOR>::transform (
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::transform (
+MappingFEField<dim,spacedim,VECTOR,DH>::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,
@@ -695,8 +658,8 @@ MappingFE<dim,spacedim,DH,VECTOR>::transform (
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
-void MappingFE<dim,spacedim,DH,VECTOR>::transform
+template<int dim, int spacedim, class VECTOR, class DH>
+void MappingFEField<dim,spacedim,VECTOR,DH>::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,
@@ -711,9 +674,9 @@ void MappingFE<dim,spacedim,DH,VECTOR>::transform
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 template < int rank >
-void MappingFE<dim,spacedim,DH,VECTOR>::transform_fields(
+void MappingFEField<dim,spacedim,VECTOR,DH>::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,
@@ -773,9 +736,9 @@ void MappingFE<dim,spacedim,DH,VECTOR>::transform_fields(
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 template < int rank >
-void MappingFE<dim,spacedim,DH,VECTOR>::transform_differential_forms(
+void MappingFEField<dim,spacedim,VECTOR,DH>::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,
@@ -807,9 +770,9 @@ void MappingFE<dim,spacedim,DH,VECTOR>::transform_differential_forms(
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 Point<spacedim>
-MappingFE<dim,spacedim,DH,VECTOR>::
+MappingFEField<dim,spacedim,VECTOR,DH>::
 transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                              const Point<dim>                                 &p) const
 {
@@ -819,6 +782,7 @@ transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_it
 
   Threads::Mutex::ScopedLock lock(mutex);
   update_internal_dofs(cell);
+
   const Quadrature<dim> point_quadrature(p);
   std::auto_ptr<InternalData>
   mdata (dynamic_cast<InternalData *> (
@@ -828,9 +792,9 @@ transform_unit_to_real_cell (const typename Triangulation<dim,spacedim>::cell_it
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 Point<spacedim>
-MappingFE<dim,spacedim,DH,VECTOR>::
+MappingFEField<dim,spacedim,VECTOR,DH>::
 transform_unit_to_real_cell_internal (const InternalData &data) const
 {
   Point<spacedim> p_real;
@@ -847,32 +811,37 @@ transform_unit_to_real_cell_internal (const InternalData &data) const
 
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 Point<dim>
-MappingFE<dim,spacedim,DH,VECTOR>::
+MappingFEField<dim,spacedim,VECTOR,DH>::
 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);
-  }
+  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;
+  try
+    {
+      initial_p_unit
+        = StaticMappingQ1<dim,spacedim>::mapping.transform_real_to_unit_cell(cell, p);
+    }
+  catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
+    {
+      // mirror the conditions of the code below to determine if we need to
+      // use an arbitrary starting point or if we just need to rethrow the
+      // exception
+      for (unsigned int d=0; d<dim; ++d)
+        initial_p_unit[d] = 0.5;
+    }
 
-
-  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);
 
+  // for (unsigned int d=0; d<dim; ++d)
+  //   initial_p_unit[d] = 0.;
+
   const Quadrature<dim> point_quadrature(initial_p_unit);
 
   UpdateFlags update_flags = update_transformation_values|update_transformation_gradients;
@@ -887,21 +856,15 @@ transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_it
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 Point<dim>
-MappingFE<dim,spacedim,DH,VECTOR>::
+MappingFEField<dim,spacedim,VECTOR,DH>::
 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();
   (void)n_shapes;
   Assert(n_shapes!=0, ExcInternalError());
@@ -911,60 +874,42 @@ transform_real_to_unit_cell_internal
   // 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
@@ -981,64 +926,47 @@ transform_real_to_unit_cell_internal
           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;
-            }
+            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;
-        }
+        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>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::compute_fill (
+MappingFEField<dim,spacedim,VECTOR,DH>::compute_fill (
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const unsigned int  n_q_points,
   const typename QProjector<dim>::DataSetDescriptor  data_set,
@@ -1047,8 +975,10 @@ MappingFE<dim,spacedim,DH,VECTOR>::compute_fill (
   std::vector<Point<spacedim> > &quadrature_points) const
 {
   const UpdateFlags update_flags(data.current_update_flags());
-  Threads::Mutex::ScopedLock lock(mutex);
-  update_internal_dofs(cell);
+  {
+    Threads::Mutex::ScopedLock lock(mutex);
+    update_internal_dofs(cell);
+  }
 
   // first compute quadrature points
   if (update_flags & update_quadrature_points)
@@ -1132,9 +1062,9 @@ MappingFE<dim,spacedim,DH,VECTOR>::compute_fill (
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::compute_fill_face (
+MappingFEField<dim,spacedim,VECTOR,DH>::compute_fill_face (
   const typename Triangulation<dim,spacedim>::cell_iterator &cell,
   const unsigned int      face_no,
   const unsigned int      subface_no,
@@ -1271,40 +1201,36 @@ MappingFE<dim,spacedim,DH,VECTOR>::compute_fill_face (
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 unsigned int
-MappingFE<dim,spacedim,DH,VECTOR>::get_degree() const
+MappingFEField<dim,spacedim,VECTOR,DH>::get_degree() const
 {
   return fe->degree;
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 ComponentMask
-MappingFE<dim,spacedim,DH,VECTOR>::get_fe_mask() const
+MappingFEField<dim,spacedim,VECTOR,DH>::get_component_mask() const
 {
   return this->fe_mask;
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 Mapping<dim,spacedim> *
-MappingFE<dim,spacedim,DH,VECTOR>::clone () const
+MappingFEField<dim,spacedim,VECTOR,DH>::clone () const
 {
-  return new MappingFE<dim,spacedim,DH,VECTOR>(*this);
+  return new MappingFEField<dim,spacedim,VECTOR,DH>(*this);
 }
 
 
-template<int dim, int spacedim, class DH, class VECTOR>
+template<int dim, int spacedim, class VECTOR, class DH>
 void
-MappingFE<dim,spacedim,DH,VECTOR>::update_internal_dofs (
+MappingFEField<dim,spacedim,VECTOR,DH>::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;
-    }
+  Assert(euler_dof_handler != 0, ExcMessage("euler_dof_handler is empty"));
 
   typename DH::cell_iterator dof_cell(*cell, euler_dof_handler);
   Assert (dof_cell->active() == true, ExcInactiveCell());
@@ -1312,27 +1238,25 @@ MappingFE<dim,spacedim,DH,VECTOR>::update_internal_dofs (
   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]);
-    }
+    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)
+template<int dim, int spacedim, class VECTOR, class DH>
+void MappingFEField<dim,spacedim,VECTOR,DH>::update_euler_vector_using_triangulation
+(VECTOR &vector) const
 {
+  AssertDimension(vector.size(), euler_dof_handler->n_dofs());
   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);
+      const Quadrature<dim> quad(fe->get_unit_support_points());
 
       MappingQ<dim,spacedim> map_q(fe->degree);
       FEValues<dim,spacedim> fe_v(map_q, *fe, quad, update_quadrature_points);
       std::vector<types::global_dof_index> dofs(fe->dofs_per_cell);
 
-      AssertDimension(fe->dofs_per_cell, support_points.size());
+      AssertDimension(fe->dofs_per_cell, fe->get_unit_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)
@@ -1343,36 +1267,45 @@ void MappingFE<dim,spacedim,DH,VECTOR>::update_euler_vector_using_triangulation
           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];
+              if (fe_mask[comp])
+                vector(dofs[q]) = points[q][fe_to_real[comp]];
             }
         }
-
     }
   else
     {
-      // Construct a MappingFE with an FEQ
+      // Construct a MappingFEField with an FEQ, and construct a
+      // standard iso-parametric interpolation on the Triangulation.
+      //
+      // Once we have this, interpolate with the given finite element
+      // to get a Mapping which is interpolatory at the support points
+      // of FE_Q(fe->degree())
       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);
+      MappingFEField<dim,spacedim,VECTOR,DH> 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();
+      const std::vector<Point<dim> > &points = feq.get_unit_support_points();
 
-      // Here construct the matrix!!!!
+      // Here construct the interpolation matrix from FE_Q^spacedim to
+      // the FiniteElement used by euler_dof_handler.
+      //
+      // The interpolation matrix is then passed to the
+      // VectorTools::interpolate() function to generate
       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]);
-            }
+          unsigned int comp_i = fe->system_to_component_index(i).first;
+          if (fe_mask[comp_i])
+            for (unsigned int j=0; j<points.size(); ++j)
+              {
+                if ( fe_to_real[comp_i] == feq.system_to_component_index(j).first)
+                  transfer(i, j) = fe->shape_value(i, points[j]);
+              }
         }
       VectorTools::interpolate(dhq, *euler_dof_handler, transfer, eulerq, vector);
     }
@@ -1382,7 +1315,7 @@ void MappingFE<dim,spacedim,DH,VECTOR>::update_euler_vector_using_triangulation
 
 
 // explicit instantiations
-#include "mapping_fe.inst"
+#include "mapping_fe_field.inst"
 
 
 DEAL_II_NAMESPACE_CLOSE
similarity index 82%
rename from source/fe/mapping_fe.inst.in
rename to source/fe/mapping_fe_field.inst.in
index fbb8b6d2839b3940c80d3787f3563309cdab261a..ce653162ef3ec53ce221c2321f40d2913d2c6f9f 100644 (file)
@@ -19,7 +19,7 @@
 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> >;
+    template class MappingFEField<deal_II_dimension, deal_II_space_dimension,  dealii::Vector<double>, dealii::DoFHandler<deal_II_dimension, deal_II_space_dimension> >;
 #endif
   }
 
diff --git a/tests/fe/mapping_fe_field_real_to_unit_b1.cc b/tests/fe/mapping_fe_field_real_to_unit_b1.cc
new file mode 100644 (file)
index 0000000..0463f9b
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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>
+#include <deal.II/fe/fe_bernstein.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_Bernstein<dim,spacedim> feb(1);
+ const FESystem<dim,spacedim> fesystem(feb, spacedim);
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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_field_real_to_unit_b2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_curved.cc
new file mode 100644 (file)
index 0000000..7efb235
--- /dev/null
@@ -0,0 +1,145 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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/fe/fe_bernstein.h>
+
+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_Bernstein<dim,spacedim> feb(2);
+ const FESystem<dim,spacedim> fesystem(feb, spacedim);
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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_field_real_to_unit_b2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.cc
new file mode 100644 (file)
index 0000000..f8e4be9
--- /dev/null
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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>
+#include <deal.II/fe/fe_bernstein.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+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 FESystem<dim,spacedim> fesystem(FE_Q<dim,spacedim>(1),1, FE_Bernstein<dim,spacedim>(2),spacedim);
+
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+
+ // Let's use only the last spacedim components
+ ComponentMask mask(spacedim+1, true);
+ mask.set(0, false);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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);
+
+      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<2,3>();
+
+  // test_real_to_unit_cell<1,3>();
+       return 0;
+}
diff --git a/tests/fe/mapping_fe_field_real_to_unit_b2_mask.output b/tests/fe/mapping_fe_field_real_to_unit_b2_mask.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
diff --git a/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b3_curved.cc
new file mode 100644 (file)
index 0000000..1ed5543
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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/fe/fe_bernstein.h>
+
+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_Bernstein<dim,spacedim> feb(3);
+ const FESystem<dim,spacedim> fesystem(feb, spacedim);
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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_field_real_to_unit_b3_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b3_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
diff --git a/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b4_curved.cc
new file mode 100644 (file)
index 0000000..2cdf1ce
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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/fe/fe_bernstein.h>
+
+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_Bernstein<dim,spacedim> feb(4);
+ const FESystem<dim,spacedim> fesystem(feb, spacedim);
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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_field_real_to_unit_b4_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b4_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
diff --git a/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_b5_curved.cc
new file mode 100644 (file)
index 0000000..7886f74
--- /dev/null
@@ -0,0 +1,146 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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/fe/fe_bernstein.h>
+
+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_Bernstein<dim,spacedim> feb(5);
+ const FESystem<dim,spacedim> fesystem(feb, spacedim);
+ DoFHandler<dim,spacedim> dhb(triangulation);
+ dhb.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhb.n_dofs());
+ const ComponentMask mask(spacedim, true);
+
+ MappingFEField<dim,spacedim> map(eulerq, dhb, 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_field_real_to_unit_b5_curved.output b/tests/fe/mapping_fe_field_real_to_unit_b5_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
similarity index 95%
rename from tests/fe/mapping_fe_real_to_unit_q1.cc
rename to tests/fe/mapping_fe_field_real_to_unit_q1.cc
index 8e0e1744a2a5f581688e18f20cb843221eafc9d1..9ac09034dc44152cf8c6d00b0294426189f74735 100644 (file)
@@ -17,7 +17,7 @@
 
 // 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
+// MappingFEField::transform_unit_to_real_cell that we get the same point as
 // we had in the beginning.
 
 #include "../tests.h"
@@ -25,7 +25,7 @@
 #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/mapping_fe_field.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/component_mask.h>
 #include <deal.II/fe/mapping_q1.h>
@@ -80,7 +80,7 @@ void test_real_to_unit_cell()
  Vector<double> eulerq(dhq.n_dofs());
  const ComponentMask mask(spacedim, true);
 
- MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+ MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
 
  map.update_euler_vector_using_triangulation(eulerq);
 
diff --git a/tests/fe/mapping_fe_field_real_to_unit_q1.output b/tests/fe/mapping_fe_field_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_field_real_to_unit_q2_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_curved.cc
new file mode 100644 (file)
index 0000000..0cce6f0
--- /dev/null
@@ -0,0 +1,147 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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(2);
+ 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);
+
+ MappingFEField<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_field_real_to_unit_q2_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q2_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
diff --git a/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.cc
new file mode 100644 (file)
index 0000000..bd8ef97
--- /dev/null
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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 FESystem<dim,spacedim> fesystem(FE_Q<dim,spacedim>(1),1, FE_Q<dim,spacedim>(2),spacedim);
+
+ DoFHandler<dim,spacedim> dhq(triangulation);
+ dhq.distribute_dofs(fesystem);
+ Vector<double> eulerq(dhq.n_dofs());
+
+ // Let's use only the last spacedim components
+ ComponentMask mask(spacedim+1, true);
+ mask.set(0, false);
+
+ MappingFEField<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_field_real_to_unit_q2_mask.output b/tests/fe/mapping_fe_field_real_to_unit_q2_mask.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_field_real_to_unit_q3_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q3_curved.cc
new file mode 100644 (file)
index 0000000..e705179
--- /dev/null
@@ -0,0 +1,148 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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(3);
+ 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);
+
+ MappingFEField<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_field_real_to_unit_q3_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q3_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
diff --git a/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc b/tests/fe/mapping_fe_field_real_to_unit_q4_curved.cc
new file mode 100644 (file)
index 0000000..3ff565d
--- /dev/null
@@ -0,0 +1,145 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2006 - 2015 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
+// MappingFEField::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_field.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>
+
+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(4);
+ 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);
+
+ MappingFEField<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_field_real_to_unit_q4_curved.output b/tests/fe/mapping_fe_field_real_to_unit_q4_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
similarity index 95%
rename from tests/fe/mapping_fe_real_to_unit_q5_curved.cc
rename to tests/fe/mapping_fe_field_real_to_unit_q5_curved.cc
index 1de8780b312dec440fad823716e65101d24889c8..1e04f9527b6f95a852bd7359166934f6b2c6fc1d 100644 (file)
@@ -17,7 +17,7 @@
 
 // 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
+// MappingFEField::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
@@ -29,7 +29,7 @@
 #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/mapping_fe_field.h>
 #include <deal.II/fe/fe_q.h>
 #include <deal.II/fe/component_mask.h>
 #include <deal.II/fe/mapping_q1.h>
@@ -37,9 +37,6 @@
 #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;
 
 
@@ -101,7 +98,7 @@ void test_real_to_unit_cell()
  Vector<double> eulerq(dhq.n_dofs());
  const ComponentMask mask(spacedim, true);
 
- MappingFE<dim,spacedim> map(eulerq, dhq, mask);
+ MappingFEField<dim,spacedim> map(eulerq, dhq, mask);
 
  map.update_euler_vector_using_triangulation(eulerq);
  
@@ -118,7 +115,6 @@ void test_real_to_unit_cell()
                  // 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());
diff --git a/tests/fe/mapping_fe_field_real_to_unit_q5_curved.output b/tests/fe/mapping_fe_field_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.