get_dof_indices on nonactive cells; FilteredIterators; Meshworker

[Scott Miller]

Hi all,

I have found that using meshworker + FilteredIterator + hanging nodes gives me the following error:

--------------------------------------------------------

An error occurred in line <3486> of file </InstalledPackages/deal.II/install.parallel/include/deal.II/dofs/dof_accessor.templates.h> in function

    void dealii::DoFCellAccessor<DH, lda>::get_dof_indices(std::vector<unsigned int>&) const [with DH = dealii::DoFHandler<2>; bool level_dof_access = false]

The violated condition was: 

    this->active()

The name and call sequence of the exception was:

    ExcMessage ("get_dof_indices() only works on active cells.")

Additional Information: 

get_dof_indices() only works on active cells.

--------------------------------------------------------

I found this in my own code, but I changed step-12 to use a FilteredIterator and see the same behavior.  I have attached the modified source for convenience.

Any ideas on this problem?  I see that there have been some recent changes to dof_accessor.templates.h (including the failed assertion above). 

Many thanks.

-Scott

[Timo Heister]

I looked at this, but I am not sure how to solve this. Basically:
meshworker/loop.h: 205
dof_info.exterior[face_no].reinit(
neighbor, neighbor->face(neighbor_face_no),
neighbor_face_no);

Here neighbor is not an active cell.
Then we end up in dof_info.h:349 and call get_indices(c) which in turn
calls get_dof_indices() even though the cell is not active.

Scott, can you test if adding
if (!c->is_level_cell() && !c->active()) return;
to dof_info.h: 296 (beginning of get_indices()) makes your code work?

Guido, why are we asking the neighbor for get_dof_indices()? Is that
needed at all?

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>



--
Timo Heister
http://www.math.tamu.edu/~heister/

[Scott Miller]

Hi Timo,

Thanks for the quick reply.  The change you suggested gets it past the current error, but 

then I get this one:

--------------------------------------------------------

An error occurred in line <1768> of file </InstalledPackages/deal.II/install.parallel/include/deal.II/lac/sparse_matrix.h> in function

    void dealii::SparseMatrix<number>::add(unsigned int, unsigned int, number) [with number = double]

The violated condition was: 

    (index != SparsityPattern::invalid_entry) || (value == 0.)

The name and call sequence of the exception was:

    ExcInvalidIndex(i, j)

Additional Information: 

The entry with index <24,22> does not exist.

--------------------------------------------------------

My guess is that get_dof_indices needs to be called to compute

the solution on the face from the neighbor.  If this is the case, then

it is definitely needed for DG elements.

-Scott

[Wolfgang Bangerth]

On 06/04/2013 03:31 PM, Scott Miller wrote:
> Hi Timo,
>
> Thanks for the quick reply. The change you suggested gets it past the
> current error, but
> then I get this one:
>
> --------------------------------------------------------
> An error occurred in line <1768> of file
> </InstalledPackages/deal.II/install.parallel/include/deal.II/lac/sparse_matrix.h>
> in function
> void dealii::SparseMatrix<number>::add(unsigned int, unsigned int,
> number) [with number = double]
> The violated condition was:
> (index != SparsityPattern::invalid_entry) || (value == 0.)
> The name and call sequence of the exception was:
> ExcInvalidIndex(i, j)
> Additional Information:
> The entry with index <24,22> does not exist.
> --------------------------------------------------------
>
> My guess is that get_dof_indices needs to be called to compute
> the solution on the face from the neighbor. If this is the case, then
> it is definitely needed for DG elements.

It sounds like it should call the function on the *child* cell, not the
inactive mother child...

W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Timo Heister]

>> My guess is that get_dof_indices needs to be called to compute
>> the solution on the face from the neighbor. If this is the case, then
>> it is definitely needed for DG elements.

Okay, it must have only "worked" with Q1 before I added that assert.

> It sounds like it should call the function on the *child* cell, not the
> inactive mother child...

Yes, which means that it was computing the wrong thing before. Glad I
removed that "feature". :-)

[Wolfgang Bangerth]

>> It sounds like it should call the function on the *child* cell, not the
>> inactive mother child...
>
> Yes, which means that it was computing the wrong thing before. Glad I
> removed that "feature". :-)

Yes :-) More testcases for the MeshWorker code would certainly help...

Best

[Scott Miller]

Just to recap:  step-12 works with "regular" iterators but not with "filtered iterators".  Should we expect the behavior to be the same? 

Why would replacing 

dof_handler.begin_active(), dof_handler.end()

with

CellFilter (IteratorFilters::Active(), dof_handler.begin_active()),

CellFilter (IteratorFilters::Active(), dof_handler.end())

cause the error we see?

-Scott

[Scott Miller]

I'm slowly getting closer to the real problem...

At the top of meshworker/loop.h, around line 32 we have:

namespace internal

{

  /**

   * Find out if an iterator supports inactive cells.

   */

  template <class DI>

  inline bool is_active_iterator(const DI &)

  {

    return false;

  }

  template <class ACCESSOR>

  inline bool is_active_iterator(const TriaActiveIterator<ACCESSOR> &)

  {

    return true;

  }

When we use a FilteredIterator, the former is_active_iterator is always called.  I think we need to change something here.  

The function call is made at precisely the location that would affect hanging nodes, around line 195.

Note that if I add another function everything seems to work just fine:

  template <class ACCESSOR>

  inline bool is_active_iterator(const FilteredIterator<TriaActiveIterator<ACCESSOR>> &)

  {

      return true;

  }

Thoughts, comments, suggestions?

-Scott

[Timo Heister]

> Note that if I add another function everything seems to work just fine:
>
> template <class ACCESSOR>
> inline bool is_active_iterator(const
> FilteredIterator<TriaActiveIterator<ACCESSOR>> &)
> {
> return true;
> }
>
> Thoughts, comments, suggestions?

That makes sense! No matter how you filter, everything will still be
active so I don't see a problem with this code.

Do you want to make a patch we can apply?

Thanks.

[Scott Miller]

Do you want to make a patch we can apply?

Sure.  It is attached.

-Scott 

[Wolfgang Bangerth]

On 06/10/2013 01:12 PM, Scott Miller wrote:
> I'm slowly getting closer to the real problem...
>
> At the top of meshworker/loop.h, around line 32 we have:
>
> namespace internal
>
> {
>
> /**
>
> * Find out if an iterator supports inactive cells.
>
> */
>
> template<classDI>
>
> inline bool is_active_iterator(const DI &)
>
> {
>
> returnfalse;
>
> }
>
>
> template <class ACCESSOR>
>
> inline bool is_active_iterator(const TriaActiveIterator<ACCESSOR> &)
>
> {
>
> returntrue;
>
> }
>
>
> When we use a FilteredIterator, the former is_active_iterator is always
> called. I think we need to change something here.

I suppose this is a question for Guido -- this function is used in one place,
namely here:

// If iterator
// is active
// and neighbor
// is refined,
// skip
// internal face.
if (internal::is_active_iterator(cell) && neighbor->has_children())
continue;

To me this suggest that what is_active_iterator(cell) should check is *whether
the cell is active*, not whether the data type of cell can only represent
active cells. However, the current implementation does the latter (it tests
the *type*, not the *object*).

Guido -- do you remember why you wrote it this way, rather than saying
cell->active()

Cheers