Error in make_hanging_node_constraints() while using parallel::distributed::Triangulation with hp::DoFHandler

[Chaitanya Dev]

Dear Deal.II Community,

I am trying to use parallel::distributed::Triangulation with hp::DoFHandler. I run into the below error in DoFTools::make_hanging_node_constraints..

I work on shape optimization using embedding domain discretization techniques. In my problems, I have to repeatedly refine the mesh and approximate the geometry. Then I assign different FESystem for the cells within and outside the geometry. 

I have implemented the code in serial with serial triangulation and hp::DoFHandler and the program works. 

Now I wish to solve larger problems so I have decided to parallelize the code and hence I use parallel::distributed::Triangulation.

Here I run into a problem in system_setup, specifically in the function DoFTools::make_hanging_node_constraints.

Hear is the error message.

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

An error occurred in line <1119> of file </calculate/temp/iwtm009/spack-stage-dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/spack-src/source/dofs/dof_tools_constraints.cc> in function

    void dealii::DoFTools::internal::make_hp_hanging_node_constraints(const DoFHandlerType&, dealii::AffineConstraints<number>&) [with DoFHandlerType = dealii::hp::DoFHandler<2, 2>; number = double]

The violated condition was: 

    cell->face(face)->child(c)->n_active_fe_indices() == 1

Additional information: 

    This exception -- which is used in many places in the library -- usually indicates that some condition which the author of the code thought must be satisfied at a certain point in an algorithm, is not fulfilled. An example would be that the first part of an algorithm sorts elements of an array in ascending order, and a second part of the algorithm later encounters an element that is not larger than the previous one.

There is usually not very much you can do if you encounter such an exception since it indicates an error in deal.II, not in your own program. Try to come up with the smallest possible program that still demonstrates the error and contact the deal.II mailing lists with it to obtain help.

Stacktrace:

-----------

#0  /opt/dealii-9.1.1/opt/spack/linux-ubuntu18.04-haswell/gcc-7.4.0/dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/lib/libdeal_II.g.so.9.1.1: void dealii::DoFTools::internal::make_hp_hanging_node_constraints<dealii::hp::DoFHandler<2, 2>, double>(dealii::hp::DoFHandler<2, 2> const&, dealii::AffineConstraints<double>&)

#1  /opt/dealii-9.1.1/opt/spack/linux-ubuntu18.04-haswell/gcc-7.4.0/dealii-9.1.1-holxtax5bwynpjaccyvwkli2ybveslrl/lib/libdeal_II.g.so.9.1.1: void dealii::DoFTools::make_hanging_node_constraints<dealii::hp::DoFHandler<2, 2>, double>(dealii::hp::DoFHandler<2, 2> const&, dealii::AffineConstraints<double>&)

#2  hp-dist-test: TestProblem<2>::setup_system()

#3  hp-dist-test: TestProblem<2>::run()

#4  hp-dist-test: main

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

This is my setup_system function.

template<int dim>

void TestProblem<dim>::setup_system() {

dof_handler.initialize(triangulation, fe_collection);

for (auto cell : dof_handler.active_cell_iterators()) {

if (!cell->is_locally_owned())

continue;

if (cell->user_flag_set())

cell->set_active_fe_index(0);

else

cell->set_active_fe_index(1);

}

dof_handler.distribute_dofs(fe_collection);

locally_owned_dofs = dof_handler.locally_owned_dofs();

DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);

hanging_node_constraints.clear();

hanging_node_constraints.reinit(locally_relevant_dofs);

DoFTools::make_hanging_node_constraints(dof_handler,

hanging_node_constraints);

hanging_node_constraints.close();

}

I did some debugging and found that the Assert fail cell->face(face)->child(c)->n_active_fe_indices() == 1 is failing on the ghost cells. To avoid this error  I tried to set_active_fe_index for the ghost cells, but this is not allowed. 

Please let me know where I am going wrong.

I have attached a MWE. You can reproduce the error with mpirun for np >= 2 cores (the code will work for np=1, because no ghost cells). 

Thank you,

Chaitanya Dev.

[David Wells]

Hi Chaitanya,

I apologize that no one has answered your question yet. Using
hp::DoFHandler with distributed triangulations is a new feature that
is not yet in a current release and there may still be some bugs in
the implementation. I suspect no one has responded to this simply
because no one knows the answer.

It's not a great solution but would it be possible to use
parallel::shared::Triangulation for your application as a stopgap
measure while we try to figure this out? That should work correctly
with hp::DoFHandler.

Thanks,
David

> --
> 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.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/626cb4ca-8971-4511-a4d1-3ed255557388%40googlegroups.com.

[Marc Fehling]

Hi Chaitanya,

This should've been fixed in #8365, which is not included in deal.II-9.1.1.

Compiling the most recent version of deal.II from the master branch should do the trick.

Marc

[Marc Fehling]

Hi Chaitanya,

we turned your minimum working example into a test case for the deal.II library #9555.

Thank you for providing your code! Would you mind to give it a look, since we reduced and changed a few parts of it?

Best,

Marc

[cd]

Hi Marc,

Thank you for providing the fix to the problem in  #8365. 

I am happy that my code was useful to setup a test case and I have no comments on the changes.

Regards,

Chaitanya Dev.