[aspect-devel] Problems with Mapping::transform_real_to_unit_cell

[Eric Heien]

Hi Wolfgang,

I've attached the parameter file that causes the problem.  It still throws the exception below, but I think this is a slightly different cause.

--------------------------------------------------------
An error occurred in line <851> of file </home/emheien/trunk/deal.II/source/grid/grid_tools.cc> in function
    std::pair<typename Container<dim, spacedim>::active_cell_iterator, dealii::Point<dim, double> > dealii::GridTools::find_active_cell_around_point(const dealii::Mapping<dim, spacedim>&, const Container<dim, spacedim>&, const dealii::Point<spacedim>&) [with int dim = 3, Container = dealii::DoFHandler, int spacedim = 3, typename Container<dim, spacedim>::active_cell_iterator = dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<3, 3> > >]
The violated condition was: 
    best_cell.first.state() == IteratorState::valid
The name and call sequence of the exception was:
    ExcPointNotFound<spacedim>(p)
Additional Information: 
The point <891669. 896973. 3.23958e+06> could not be found inside any of the subcells of a coarse grid cell.

Stacktrace:
-----------
#0  /home/emheien/trunk/deal.II/lib/libdeal_II.g.so.7.2.pre: std::pair<dealii::DoFHandler<3, 3>::active_cell_iterator, dealii::Point<3, double> > dealii::GridTools::find_active_cell_around_point<3, dealii::DoFHandler, 3>(dealii::Mapping<3, 3> const&, dealii::DoFHandler<3, 3> const&, dealii::Point<3, double> const&)
#1  /home/emheien/trunk/deal.II/lib/libdeal_II.g.so.7.2.pre: dealii::Functions::FEFieldFunction<3, dealii::DoFHandler<3, 3>, dealii::TrilinosWrappers::MPI::BlockVector>::compute_point_locations(std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > > const&, std::vector<dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<3, 3> > >, std::allocator<dealii::TriaActiveIterator<dealii::DoFCellAccessor<dealii::DoFHandler<3, 3> > > > >&, std::vector<std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >, std::allocator<std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > > > >&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&) const
#2  /home/emheien/trunk/deal.II/lib/libdeal_II.g.so.7.2.pre: dealii::Functions::FEFieldFunction<3, dealii::DoFHandler<3, 3>, dealii::TrilinosWrappers::MPI::BlockVector>::vector_value_list(std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > > const&, std::vector<dealii::Vector<double>, std::allocator<dealii::Vector<double> > >&) const
#3  ./lib/aspect: aspect::Postprocess::ParticleSet<3>::get_particle_velocities(dealii::parallel::distributed::Triangulation<3, 3> const&, dealii::TrilinosWrappers::MPI::BlockVector const&, dealii::DoFHandler<3, 3> const&, dealii::Mapping<3, 3> const&, std::vector<dealii::Point<3, double>, std::allocator<dealii::Point<3, double> > >&)
#4  ./lib/aspect: aspect::Postprocess::ParticleSet<3>::advect_particles(dealii::parallel::distributed::Triangulation<3, 3> const&, double, dealii::TrilinosWrappers::MPI::BlockVector const&, dealii::DoFHandler<3, 3> const&, dealii::Mapping<3, 3> const&)
#5  ./lib/aspect: aspect::Postprocess::ParticleSet<3>::execute(dealii::TableHandler&)
#6  ./lib/aspect: aspect::Postprocess::Manager<3>::execute(dealii::TableHandler&)
#7  ./lib/aspect: aspect::Simulator<3>::postprocess()
#8  ./lib/aspect: aspect::Simulator<3>::run()
#9  ./lib/aspect: main


-Eric

On Jun 27, 2012, at 12:00 PM, Wolfgang Bangerth wrote:

> 
> Folks,
> coming back to this: I think I've solved this now, hopefully correctly: 
> whenever the iteration in Mapping::transform_real_to_unit_cell fails, 
> then we throw an exception. I've added a line search to the iteration, 
> so I think it is reasonable to assume that if the iteration fails, then 
> because the point is, in fact, outside the cell, and I've modified the 
> various callers of the function to expect that an exception may be 
> thrown. Thanks for all of your help in figuring out where the problem 
> actually lies!
> 
> Eric: Do you still have the input file that triggered the problem? If 
> so, can you (i) send it to me, and (ii) verify that the problem is 
> indeed gone now?
> 
> Best
>  W.
> 
> ------------------------------------------------------------------------
> Wolfgang Bangerth               email:            bangerth at math.tamu.edu
>                                 www: http://www.math.tamu.edu/~bangerth/
> 
> 
> 
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel

-------------- next part --------------
A non-text attachment was scrubbed...
Name: particle_q4_test.prm
Type: application/octet-stream
Size: 13738 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/aspect-devel/attachments/20120629/fe1914e3/attachment.obj