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