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
(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@>.
/*@{*/
/**
- * 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);
* 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;
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,
* 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
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;
};
* 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
* 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;
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,
/**
- * 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;
};
/*@}*/
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,
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,
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(),
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(),
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(),
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(),
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;
}
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
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
#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>
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),
-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()),
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)
{
}
-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),
{}
-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));
-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.
// 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,
-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();
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);
}
-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,
-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));
}
-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));
}
-// 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,
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;
-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,
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);
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,
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,
-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,
}
-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,
}
-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,
}
-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,
-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
{
Threads::Mutex::ScopedLock lock(mutex);
update_internal_dofs(cell);
+
const Quadrature<dim> point_quadrature(p);
std::auto_ptr<InternalData>
mdata (dynamic_cast<InternalData *> (
}
-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;
-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;
}
-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());
// 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
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,
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)
}
-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,
}
-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());
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)
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);
}
// explicit instantiations
-#include "mapping_fe.inst"
+#include "mapping_fe_field.inst"
DEAL_II_NAMESPACE_CLOSE
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
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
// 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"
#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>
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);
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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
// 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
#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>
#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;
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);
// 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());
--- /dev/null
+
+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