* equals the given point @p p. If this is the case then this function
* throws an exception of type Mapping::ExcTransformationFailed . Whether
* the given point @p p lies outside the cell can therefore be determined by
- * checking whether the return reference coordinates lie inside or outside
+ * checking whether the returned reference coordinates lie inside or outside
* the reference cell (e.g., using GeometryInfo::is_inside_unit_cell()) or
* whether the exception mentioned above has been thrown.
*
*/
virtual
std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
- get_vertices
- (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+ get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
/**
* Return a pointer to a copy of the present object. The caller of this copy
const VECTOR &euler_vector,
const DoFHandler<dim,spacedim> &euler_dof_handler) DEAL_II_DEPRECATED;
+ /**
+ * Return the mapped vertices of the cell. For the current class, this function does
+ * not use the support points from the geometry of the current cell but
+ * instead evaluates an externally given displacement field in addition to
+ * the geometry of the cell.
+ */
+ virtual
+ std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+ get_vertices (const typename Triangulation<dim,spacedim>::cell_iterator &cell) const;
+
/**
* Return a pointer to a copy of the present object. The caller of this copy
* then assumes ownership of it.
transform_real_to_unit_cell (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
const Point<spacedim> &p) const
{
- // first a Newton iteration based on a Q1 mapping to get a good starting
- // point, the idea being that this is cheaper than trying to start with the
- // real mapping and likely also more robust.
- //
- // that said, this doesn't always work: there are cases where the point is
- // outside the cell and the inverse mapping doesn't converge. in that case,
- // use the center point of the cell as a starting point if we are to go on
- // using the full mapping, or just propagate up the exception if we had no
- // intention of continuing with the full mapping
- Point<dim> initial_p_unit;
- try
- {
- initial_p_unit = q1_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
- if (cell->has_boundary_lines()
- ||
- use_mapping_q_on_all_cells
- ||
- (dim!=spacedim) )
- {
- for (unsigned int d=0; d<dim; ++d)
- initial_p_unit[d] = 0.5;
- }
- else
- throw;
- }
-
- // then a Newton iteration based on the full MappingQ if we need this. note
- // that for interior cells with dim==spacedim, the mapping used is in fact a
- // Q1 mapping, so there is nothing we need to do unless the iteration above
- // failed
if (cell->has_boundary_lines()
||
use_mapping_q_on_all_cells
||
(dim!=spacedim) )
- {
- // 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);
-
- return this->transform_real_to_unit_cell_internal(cell, p, initial_p_unit);
- }
+ return MappingQGeneric<dim,spacedim>::transform_real_to_unit_cell(cell, p);
else
- return initial_p_unit;
+ return q1_mapping->transform_real_to_unit_cell(cell, p);
}
// .... COMPUTE MAPPING SUPPORT POINTS
+template <int dim, class EulerVectorType, int spacedim>
+std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell>
+MappingQEulerian<dim, EulerVectorType, spacedim>::
+get_vertices
+(const typename Triangulation<dim,spacedim>::cell_iterator &cell) const
+{
+ // get the vertices as the first 2^dim mapping support points
+ std::vector<Point<spacedim> > a;
+ compute_mapping_support_points(cell, a);
+
+ std_cxx11::array<Point<spacedim>, GeometryInfo<dim>::vertices_per_cell> vertex_locations;
+ std::copy (a.begin(),
+ a.begin()+GeometryInfo<dim>::vertices_per_cell,
+ vertex_locations.begin());
+
+ return vertex_locations;
+}
+
+
+
template <int dim, class EulerVectorType, int spacedim>
void
MappingQEulerian<dim, EulerVectorType, spacedim>::
compute_mapping_support_points (const typename Triangulation<dim,spacedim>::cell_iterator &cell,
std::vector<Point<spacedim> > &a) const
{
-
// first, basic assertion with respect to vector size,
const types::global_dof_index n_dofs = euler_dof_handler->n_dofs();
#include <deal.II/fe/fe.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/mapping_q1.h>
#include <cmath>
#include <algorithm>
}
- // Find the initial value for the Newton iteration by a normal
- // projection to the least square plane determined by the vertices
- // of the cell
- std::vector<Point<spacedim> > a;
- compute_mapping_support_points (cell,a);
- Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
- ExcInternalError());
- const Point<dim> initial_p_unit =
- internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+ Point<dim> initial_p_unit;
+ if (polynomial_degree == 1)
+ {
+ // Find the initial value for the Newton iteration by a normal
+ // projection to the least square plane determined by the vertices
+ // of the cell
+ std::vector<Point<spacedim> > a;
+ compute_mapping_support_points (cell,a);
+ Assert(a.size() == GeometryInfo<dim>::vertices_per_cell,
+ ExcInternalError());
+ initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+ }
+ else
+ {
+ try
+ {
+ // Find the initial value for the Newton iteration by a normal
+ // projection to the least square plane determined by the vertices
+ // of the cell
+ //
+ // we do this by first getting all support points, then
+ // throwing away all but the vertices, and finally calling
+ // the same function as above
+ std::vector<Point<spacedim> > a;
+ compute_mapping_support_points (cell,a);
+ a.resize(GeometryInfo<dim>::vertices_per_cell);
+ initial_p_unit = internal::MappingQ1::transform_real_to_unit_cell_initial_guess<dim,spacedim>(a,p);
+ }
+ catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ initial_p_unit[d] = 0.5;
+ }
+
+ // 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);
+ }
// perform the Newton iteration and return the result. note that
// this statement may throw an exception, which we simply pass up to