get_dof_indices() only works on active cells - blocks vector

[Juan Felipe Giraldo]

Hi all!

I am currently working on the linear advection problem by using a DG formulation, similar to the problem in step 12 https://www.dealii.org/current/doxygen/deal.II/step_12.html, but using a mixed formulation as it is shown in step 20 https://www.dealii.org/current/doxygen/deal.II/step_20.html. I can solve the problem without problems using a global refinement, however, when I implement an adaptative refinement I have the following problem after the 1st refinement cycle

An error occurred in line <3792> of file </home/arsen/Documents/instaladores/dealii/include/deal.II/dofs/dof_accessor.templates.h> in function

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

The violated condition was: 

    this->active()

Additional information: 

    get_dof_indices() only works on active cells.

 

I have seen that some people had this problem a couple of years ago, but I have not seen a proper solution to solve it. 

Anyone knows how can I solve it?

Thank you so much,

Regards,

Juan Felipe Giraldo

[David Wells]

Dear Juan,

My best guess is that, at some point in the loop over all active
cells, you are also doing some computations involving the neighbor of
each cell. TriaAccessor::neighbor() returns the neighboring cell which
has, at most, the same level number as the current cell: i.e., if the
neighbor cell has been refined but the current cell has not, then
trying to get the degrees of freedom from the neighbor will throw
exactly this exception since the cell returned by the call to
neighbor() is inactive.

This is only a guess. Could you show us the part of your code that
triggers this exception? That will help us give you a better answer :)

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/38d23dcb-3028-45eb-a539-737e72c597b0%40googlegroups.com.

[Juan Felipe Giraldo]

Dear David,

 

Thank you for your quick reply. 

Yes, I need the neighbour information of each cell to compute the boundary integration in the DG formulation. Therefore, I get the neighbouring cell in the loop of the active cell by using "const auto neighbor = cell->neighbor(face_n)";  (similar to step 30). 

I am quite sure that the problem is what you are describing, thus it occurs only when I use the adaptative refinement. 

The problem appears here:

#0  ./FEM-DGdeal: dealii::DoFCellAccessor<dealii::DoFHandler<2, 2>, false>::get_dof_indices(std::vector<unsigned int, std::allocator<unsigned int> >&) const

#1  ./FEM-DGdeal: FEMwDG<2>::assemble_system(unsigned int)

#2  ./FEM-DGdeal: FEMwDG<2>::run(unsigned int, unsigned int, unsigned int, unsigned int)

#3  ./FEM-DGdeal: main

And the part of the code where I call the get_dof_indices is this:

    void FEMwDG<dim>::assemble_system(const unsigned int typecase)

    {

      const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;

      std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

      std::vector<types::global_dof_index> neighbor_dof_indices(dofs_per_cell);

      const UpdateFlags update_flags = update_values | update_gradients |

                                       update_quadrature_points |

                                       update_JxW_values;

      const UpdateFlags face_update_flags =  update_values | update_quadrature_points | update_JxW_values |

                                             update_normal_vectors;

      const UpdateFlags neighbor_face_update_flags = update_values;

      FEValues<dim>     fe_values(mapping, fe, quadrature, update_flags);

      FEFaceValues<dim> fe_face_values(mapping, fe, face_quadrature, face_update_flags);

      FESubfaceValues<dim> fe_subface_values(mapping,  fe, face_quadrature, face_update_flags);

      FEFaceValues<dim>  fe_face_neighbor(mapping,  fe, face_quadrature,neighbor_face_update_flags);

      FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell); //matrix dominio

      FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);

      FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);

      FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);

      Vector<double>     cell_vector(dofs_per_cell);

      const FEValuesExtractors::Scalar verror(0);

      const FEValuesExtractors::Scalar variable(1);

      for (const auto &cell : dof_handler.active_cell_iterators())

        {

          fe_values.reinit(cell);   ui_vi_matrix = 0;  cell_vector    = 0;

      //--------------------------------------DOMAIN INTEGRATOR--------------------------------------------------------------------------

             adv.assemble_cell_term(fe_values, ui_vi_matrix, cell_vector, verror, variable,cell, typecase);

             cell->get_dof_indices(local_dof_indices);

                for (unsigned int face_n = 0;  face_n < GeometryInfo<dim>::faces_per_cell;  ++face_n)

                  {

                    //-------------------------------------FACE AT THE EXTERNAL BOUNDARY-----------------------------------------------------------------------------

                    const auto face = cell->face(face_n);

                     if (face->at_boundary())

                       {

                         fe_face_values.reinit(cell, face_n);

                         adv.assemble_boundary_term(fe_face_values, ui_vi_matrix, cell_vector, verror, variable);

                       }

                    //-------------------------------------FACE AT THE INTERNAL BOUNDARY-----------------------------------------------------------------------------

                   else

                   {

                      const auto neighbor = cell->neighbor(face_n);

                         if ((dim > 1) && (cell->index() < neighbor->index()))

                      {

                         const unsigned int neighbor2 =cell->neighbor_face_no(face_n);

                         ue_vi_matrix = 0; ui_ve_matrix = 0;  ue_ve_matrix = 0;

                         fe_face_values.reinit(cell, face_n);

                         fe_face_neighbor.reinit(neighbor, neighbor2);

                         adv.assemble_face_term(fe_face_values, fe_face_neighbor,  ui_vi_matrix, ue_vi_matrix, ui_ve_matrix, ue_ve_matrix,  verror,  variable,  typecase);

                         neighbor->get_dof_indices(neighbor_dof_indices);

                     // ---------------LOCAL ASSEMBLE ----------------------------------------------------------

                        for (unsigned int i = 0; i < dofs_per_cell; ++i)

                           for (unsigned int j = 0; j < dofs_per_cell; ++j)

                             {

                                 system_matrix.add(local_dof_indices[i],  neighbor_dof_indices[j],  ue_vi_matrix(i, j));

                                 system_matrix.add(neighbor_dof_indices[i],  local_dof_indices[j],  ui_ve_matrix(i, j));

                                 system_matrix.add(neighbor_dof_indices[i], neighbor_dof_indices[j], ue_ve_matrix(i, j));

                           }

                      }

                   }

             }

        // ----------------------------------GLOBAL ASSEMBLE ----------------------------------------------------------

                   for (unsigned int i = 0; i < dofs_per_cell; ++i)

                     {for (unsigned int j = 0; j < dofs_per_cell; ++j)

                       {system_matrix.add(local_dof_indices[i],

                                          local_dof_indices[j],

                                          ui_vi_matrix(i, j));}

                        system_rhs(local_dof_indices[i]) += cell_vector(i);}

        }

  }

Thanks :),

Juan Giraldo

> To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.

[luca.heltai]

Dear Juan,

The logic you should use is: when your neighbour is finer, you should not assemble the full face term from the coarse neighbour. Instead, you have to assemble your DG terms using a FESubfaceValues initialized with the correct arguments on the coarse face, to match *each* of the fine FEFaceValues. This logic is tricky to get right.

May I suggest that you use the MeshWorker::mesh_loop function?

That function takes care of these issues for you, by calling a user provided assembly function for cells, faces and boundary faces, and passing to that function the right arguments with which you can initialise your FEFaceValues and FESubFaceValues.

You can see it in action in step-12 (development version of deal.II), or in the tests/meshworker/mesh_loop_*.cc files and tests/meshworker/scratch_data_*.cc

Best,
Luca.

> 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/3d98145f-d1ed-4b4d-ae62-22d61df5ebe9%40googlegroups.com.

[Juan Felipe Giraldo]

Dear Luca,

Thank you very much for your reply, your answer was very precise and, therefore, I could solve my problem. Effectively, I was not considering the subfaces when I was assembling, and of course the "get_dof_indices" function was not recognizing the nodes neighbour nodes those cells.

Note: As I am working with mixed formulation, I need to compute the cell, boundary and face terms for two different variables (in blocks). I had some problems to declare them by separate using the meshworker. So, it was easier for me to declare the functions independently.

Regards,

Juan Giraldo 

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/3d98145f-d1ed-4b4d-ae62-22d61df5ebe9%40googlegroups.com.