Re: [deal.II] Bug in GridTools::find_active_cell_around_point ?

Steffen Brinckmann Fri, 25 May 2012 00:15:47 -0700

Dear Timo,

in the mean-time I was able to figure out the problem, but not able topost to the deal.ii list.


The problem is GridTools::find_active_cell_around_point:
1. it find the closest vertex

2. and then loops through all cells to identify the correct cell that isconnected to this vertex

This method works fine for all uniformly refined grids.

However, if you have a locally refined grid with hanging nodes, than theclosest vertex might be a hanging node, and your point is in the"hanging-element". The loop through the cells will fail, because thepoint is in a cell, which has this hanging node not as vertex.

I replaced the GridTools::find_active_cell_around_point with a simpleloop through all cells to identify the correct one using point_inside().The hierachical search method of going from the parents to the childrenalso works fine.

I can certainly try to reduce the code to prove the issue, but I guessthe above argument might be sufficient.

Steffen


On 05/24/2012 10:03 PM, Timo Heister wrote:

Hi Steffen,

can you please try to reduce this to a minimal test case that we can
study (just by generating the grid and calling point_inside() with the
problematic coordinates)?

Thanks,
Timo

On Fri, May 11, 2012 at 1:12 PM, Steffen Brinckmann
<[email protected]>  wrote:

Dear All,
I am not able to directly localize the - what I would call - bug, but here
is the description on how far I got.
1. I test whether a point is inside my domain with a loop across all cells
and the "cell->point_inside (point)" command. Yes my point is inside, also
visual inspection and geometry: the point is inside.

2. I run  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&)
via the FEFieldFunction and that raises the exception
------------------------------------------------------
The violated condition was:
    best_cell.first.state() == IteratorState::valid
The name and call sequence of the exception was:
    ExcPointNotFound<dim>(p)
Additional Information:
The point<304.767 -57.0113 254.766>  could not be found inside any of the
subcells of a coarse grid cell.
-------------------------------------------------------

3. I rerun, and the exception is repeatable: every time I run it.

Running a debugger I found that the exception is thrown after (not during)
the return statement of
-------------------------------------------------------
template<typename T, typename P>
inline
T&  SmartPointer<T,P>::operator * () const
{
  Assert(t != 0, ExcNotInitialized());
  return *t;
}
-------------------------------------------------------

Can someone help?
THANKS,
Steffen



--------------------------------------------------------
An error occurred in line<830>  of file<...
/deal.II/source/grid/grid_tools.cc>  in function
    std::pair<typename Container<dim, spacedim>::active_cell_iterator,
dealii::Point<spacedim>  >
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<dim>(p)
Additional Information:
The point<304.767 -57.0113 254.766>  could not be found inside any of the
subcells of a coarse grid cell.

Stacktrace:
-----------
#0   deal.II/lib/libdeal_II.g.so.7.1.0: 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  deal.II/lib/libdeal_II.g.so.7.1.0: dealii::Functions::FEFieldFunction<3,
dealii::DoFHandler<3, 3>, dealii::Vector<double>

::vector_gradient(dealii::Point<3, double>  const&,

std::vector<dealii::Tensor<1, 3, double>, std::allocator<dealii::Tensor<1,
3, double>  >  >&) const
#2  exe: ares::FEM::getStress(dealii::Point<3, double>)
#3  exe: ares::FEM::getStress(Eigen::Matrix<float, 3, 1, 0, 3, 1>)

_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii



--
Dr. Steffen Brinckmann
Groupleader
ICAMS
Ruhr-University Bochum
Stiepeler Str. 129
44780 Bochum
Germany

Tel.   +49 234 32 29372
Fax    +49 234 32 14977
Skype: sbrinckm
Email: [email protected]

_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii

Re: [deal.II] Bug in GridTools::find_active_cell_around_point ?

Reply via email to