<ol>
+<li> Fixed: The MappingQ1::transform_real_to_unit_cell function as
+well as the equivalent ones in derived classes sometimes get into
+trouble if they are asked to compute the preimage of this point
+in reference cell coordinates. This is because for points outside
+the reference cell, the mapping from unit to real cell is not
+necessarily invertible, and consequently the Newton iteration to
+find the preimage did not always converge, leading to an exception.
+While this is not entirely wrong (we could, after all, not compute
+the desired quantity), not all callers of this function were prepared
+to accept this result -- in particular, the function
+CellAccessor<3>::point_inside should have really just returned false
+in such cases but instead let the exception so generated propagate
+through. This should now be fixed.
+<br>
+(Wolfgang Bangerth, Eric Heien, Sebastian Pauletti, 2012/06/27)
+
<li> Fixed: The function VectorTools::compute_no_normal_flux_constraints had
a bug that led to an exception whenever we were computing constraints for
vector fields located on edges shared between two faces of a 3d cell if those
/**
* Transforms the point @p p on
- * the real cell to the point
- * @p p_unit on the unit cell
- * @p cell and returns @p p_unit.
+ * the real @p cell to the corresponding
+ * point on the unit cell, and
+ * return its coordinates.
*
* In the codimension one case,
* this function returns the
* point @p p on the curve or
* surface identified by the @p
* cell.
+ *
+ * @note Polynomial mappings from
+ * the reference (unit) cell coordinates
+ * to the coordinate system of a real
+ * cell are not always invertible if
+ * the point for which the inverse
+ * mapping is to be computed lies
+ * outside the cell's boundaries.
+ * In such cases, the current function
+ * may fail to compute a point on
+ * the reference cell whose image
+ * under the mapping 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 of outside the reference
+ * cell (e.g., using
+ * GeometryInfo::is_inside_unit_cell)
+ * or whether the exception mentioned
+ * above has been thrown.
*/
virtual Point<dim>
transform_real_to_unit_cell (
*/
DeclException0 (ExcInvalidData);
+
+ /**
+ * Computing the mapping between a
+ * real space point and a point
+ * in reference space failed, typically because the given point
+ * lies outside the cell where the inverse mapping is not
+ * unique.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException0(ExcTransformationFailed);
+
private:
/**
* point @p p on the curve or
* surface identified by the @p
* cell.
+ *
+ * @note Polynomial mappings from
+ * the reference (unit) cell coordinates
+ * to the coordinate system of a real
+ * cell are not always invertible if
+ * the point for which the inverse
+ * mapping is to be computed lies
+ * outside the cell's boundaries.
+ * In such cases, the current function
+ * may fail to compute a point on
+ * the reference cell whose image
+ * under the mapping 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 of outside the reference
+ * cell (e.g., using
+ * GeometryInfo::is_inside_unit_cell)
+ * or whether the exception mentioned
+ * above has been thrown.
*/
virtual Point<dim>
transform_real_to_unit_cell (
* point @p p on the curve or
* surface identified by the @p
* cell.
+ *
+ * @note Polynomial mappings from
+ * the reference (unit) cell coordinates
+ * to the coordinate system of a real
+ * cell are not always invertible if
+ * the point for which the inverse
+ * mapping is to be computed lies
+ * outside the cell's boundaries.
+ * In such cases, the current function
+ * may fail to compute a point on
+ * the reference cell whose image
+ * under the mapping 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 of outside the reference
+ * cell (e.g., using
+ * GeometryInfo::is_inside_unit_cell)
+ * or whether the exception mentioned
+ * above has been thrown.
*/
virtual Point<dim>
transform_real_to_unit_cell (
virtual
bool preserves_vertex_locations () const;
- /**
- * The Newton iteration computing
- * the mapping between a
- * real space point and a point
- * in reference space did not
- * converge, thus resulting in a
- * wrong mapping.
- *
- * @ingroup Exceptions
- */
- DeclException0(ExcTransformationFailed);
-
-
protected:
/* Trick to templatize transform_real_to_unit_cell<dim, dim+1> */
template<int dim_>
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.
- Point<dim> initial_p_unit =
- MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
+ // 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
+ Point<dim> initial_p_unit;
+ bool initial_p_unit_is_valid = false;
+ try
+ {
+ initial_p_unit
+ = MappingQ1<dim,spacedim>::transform_real_to_unit_cell(cell, p);
+
+ initial_p_unit_is_valid = true;
+ }
+ catch (const typename Mapping<dim,spacedim>::ExcTransformationFailed &)
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ initial_p_unit[d] = 0.5;
+ }
// 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
- if (cell->has_boundary_lines() ||
- use_mapping_q_on_all_cells ||
+ // do unless the iteration above failed
+ if ((initial_p_unit_is_valid == false)
+ ||
+ cell->has_boundary_lines()
+ ||
+ use_mapping_q_on_all_cells
+ ||
(dim!=spacedim) )
{
// use the full mapping. in case the
DEAL_II_NAMESPACE_OPEN
-
-
-
template <int dim, int spacedim>
const unsigned int MappingQ1<dim,spacedim>::n_shape_functions;
// Ignore non vertex support points.
mdata->mapping_support_points.resize(GeometryInfo<dim>::vertices_per_cell);
- // perform the Newton iteration and
- // return the result
+ // perform the Newton iteration and
+ // return the result. note that this
+ // statement may throw an exception, which
+ // we simply pass up to the caller
return transform_real_to_unit_cell_internal(cell, p, initial_p_unit,
*mdata);
}
Point<spacedim> f = p_real-p;
const double eps = 1.e-12*cell->diameter();
- const unsigned int loop_limit = 50;
+ const unsigned int newton_iteration_limit = 20;
- unsigned int loop=0;
- while (f.square()>eps*eps && loop++<loop_limit)
+ unsigned int newton_iteration=0;
+ while (f.square()>eps*eps)
{
// f'(x)
Tensor<2,spacedim> df;
}
// Solve [f'(x)]d=f(x)
- Tensor<1,spacedim> d;
- contract (d, invert(df), static_cast<const Tensor<1,spacedim>&>(f));
-
- // update of p_unit. The
- // spacedimth component of
- // transformed point is simply
- // ignored in codimension one
- // case. When this component is
- // not zero, then we are
- // projecting the point to the
- // surface or curve identified
- // by the cell.
- for(unsigned int i=0;i<dim; ++i)
- p_unit[i] -= d[i];
-
- // shape values and derivatives
- // at new p_unit point
- compute_shapes(std::vector<Point<dim> > (1, p_unit), mdata);
-
- // f(x)
- p_real = transform_unit_to_real_cell_internal(mdata);
- f = p_real-p;
- }
+ Tensor<1,spacedim> delta;
+ contract (delta, invert(df), static_cast<const Tensor<1,spacedim>&>(f));
- // Here we check that in the last
- // execution of while the first
- // condition was already wrong,
- // meaning the residual was below
- // eps. Only if the first condition
- // failed, loop will have been
- // increased and tested, and thus
- // havereached the limit.
- AssertThrow (loop<loop_limit, ExcTransformationFailed());
+ // do a line search
+ double step_length = 1;
+ do
+ {
+ // update of p_unit. The
+ // spacedimth component of
+ // transformed point is simply
+ // ignored in codimension one
+ // case. When this component is
+ // not zero, then we are
+ // projecting the point to the
+ // surface or curve identified
+ // by the cell.
+ Point<dim> p_unit_trial = p_unit;
+ for(unsigned int i=0;i<dim; ++i)
+ p_unit_trial[i] -= step_length * delta[i];
+
+ // shape values and derivatives
+ // at new p_unit point
+ compute_shapes(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_real_trial-p;
+
+ // 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() < f.norm())
+ {
+ p_real = p_real_trial;
+ p_unit = p_unit_trial;
+ f = f_trial;
+ break;
+ }
+ else if (step_length > 0.05)
+ step_length /= 2;
+ else
+ goto failure;
+ }
+ while (true);
+
+ ++newton_iteration;
+ if (newton_iteration > newton_iteration_limit)
+ 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>();
}
}
}
+//TODO: implement a line search here in much the same way as for
+// the corresponding function above that does so for codim==0.
p_minus_F = p;
p_minus_F -= transform_unit_to_real_cell_internal(mdata);
// failed, loop will have been
// increased and tested, and thus
// have reached the limit.
- AssertThrow (loop<loop_limit, ExcTransformationFailed());
+ AssertThrow (loop<loop_limit, (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
return p_unit;
}
if ((p[d] < minp[d]) || (p[d] > maxp[d]))
return false;
- // now we need to check more
- // carefully: transform to the
- // unit cube
- // and check there.
- const TriaRawIterator< CellAccessor<dim,spacedim> > cell_iterator (*this);
- return (GeometryInfo<dim>::is_inside_unit_cell (
- StaticMappingQ1<dim,spacedim>::mapping.transform_real_to_unit_cell(cell_iterator, p)));
+ // now we need to check more carefully: transform to the
+ // unit cube and check there. unfortunately, this isn't
+ // completely trivial since the transform_real_to_unit_cell
+ // function may throw an exception that indicates that the
+ // point given could not be inverted. we take this as a sign
+ // that the point actually lies outside, as also documented
+ // for that function
+ try
+ {
+ const TriaRawIterator< CellAccessor<dim,spacedim> > cell_iterator (*this);
+ return (GeometryInfo<dim>::is_inside_unit_cell
+ (StaticMappingQ1<dim,spacedim>::mapping.transform_real_to_unit_cell(cell_iterator, p)));
+ }
+ catch (const Mapping<dim,spacedim>::ExcTransformationFailed &e)
+ {
+ return false;
+ }
}
MappingQ<dim> map(4);
typename Triangulation<dim >::active_cell_iterator
cell = triangulation.begin_active();
- deallog << map.transform_real_to_unit_cell(cell,p) << std::endl;
+
+ try
+ {
+ map.transform_real_to_unit_cell(cell,p);
+ }
+ catch (const typename Mapping<dim>::ExcTransformationFailed &)
+ {
+ deallog << "Point is outside!" << std::endl;
+ }
}
-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::Point is outside!
MappingQ1<dim> map;
typename Triangulation<dim >::active_cell_iterator
cell = triangulation.begin_active();
- deallog << map.transform_real_to_unit_cell(cell,p) << std::endl;
+ try
+ {
+ map.transform_real_to_unit_cell(cell,p);
+ }
+ catch (const typename Mapping<dim>::ExcTransformationFailed &)
+ {
+ deallog << "Point is outside!" << std::endl;
+ }
}
-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::Point is outside!
+++ /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
//-----------------------------------------------------------------------------
-// a redux of mapping_real_to_unit_q4_sphere_x.
-// It seems that the point is outside the cell??
-
+// a redux of mapping_real_to_unit_q4_sphere. one doesn't in fact need
+// the Q4 mapping, the problem already happens in the initial guess
+// generation using a Q1 mapping
#include "../tests.h"
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/fe/mapping_q.h>
-#include <deal.II/numerics/data_out.h>
void test_real_to_unit_cell()
{
triangulation.create_triangulation (vertices, cells,
SubCellData());
-
- FE_Q<dim> fe (1);
- DoFHandler<dim> dh(triangulation);
- dh.distribute_dofs (fe);
- Vector<double> dummy;
- dummy.reinit(dh.n_dofs());
- DataOut<dim, DoFHandler<dim> > data_out;
- data_out.attach_dof_handler (dh);
- data_out.add_data_vector (dummy,"dummy");
- data_out.build_patches ();
- std::ofstream output ("mapping_real_to_unit_q4_sphere_xx/plot.vtk");
- data_out.write_vtk (output);
-
-
- // const Point<dim> p (-3.56413e+06, 1.74215e+06, 2.14762e+06);
- // MappingQ1<dim> map;
- // typename Triangulation<dim >::active_cell_iterator
- // cell = triangulation.begin_active();
- // deallog << map.transform_real_to_unit_cell(cell,p) << std::endl;
+ // set the boundary indicator for
+ // one face and adjacent edges of
+ // the single cell
+ triangulation.set_boundary (1, boundary);
+ triangulation.begin_active()->face(5)->set_all_boundary_indicators (1);
+
+ // now try to find the coordinates
+ // of the following point in the
+ // reference coordinate system of
+ // the cell
+ const Point<dim> p (-3.56413e+06, 1.74215e+06, 2.14762e+06);
+ MappingQ1<dim> map;
+ typename Triangulation<dim >::active_cell_iterator
+ cell = triangulation.begin_active();
+
+ // the following call will fail
+ // because the point is outside the
+ // bounds of the cell. catch the
+ // exception and make sure that
+ // there is output that documents
+ // that the function threw an
+ // exception
+ try
+ {
+ map.transform_real_to_unit_cell(cell,p);
+ }
+ catch (const std::exception &e)
+ {
+ deallog << "Caught exception of type " << typeid(e).name()
+ << std::endl;
+ }
}
int
main()
{
- std::ofstream logfile ("mapping_real_to_unit_q4_sphere_xx/output");
+ std::ofstream logfile ("mapping_real_to_unit_q4_sphere_y/output");
deallog.attach(logfile);
deallog.depth_console(0);
deallog.threshold_double(1.e-10);
--- /dev/null
+
+DEAL::Caught exception of type N6dealii7MappingILi3ELi3EE23ExcTransformationFailedE
projects / dealii-svn.git / commitdiff