Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Error when extract face dofs of face boardering FE_Nothing #13418

Closed
simonsticko opened this issue Feb 18, 2022 · 2 comments · Fixed by #13430
Closed

Error when extract face dofs of face boardering FE_Nothing #13418

simonsticko opened this issue Feb 18, 2022 · 2 comments · Fixed by #13430

Comments

@simonsticko
Copy link
Contributor

I thought the following would work, but it results in an exception. Is this a bug or not?

Set up a triangulation with 2 cells:

|cell 0|cell 1|

and a DoFHandler with elements over the triangulation:

|FE_Q|FE_Nothing|

Try to extract the face dofs of cell 0 at the shared face.

#include <deal.II/dofs/dof_handler.h>

#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>

#include <deal.II/grid/grid_generator.h>

#include <vector>

using namespace dealii;

template <int dim>
void
make_two_element_mesh(Triangulation<dim> &triangulation)
{
  const Point<dim> lower_left;
  Point<dim>       upper_right;

  const unsigned int        n_cells = 2;
  std::vector<unsigned int> repetitions;
  upper_right[0] = n_cells;
  repetitions.push_back(n_cells);
  for (unsigned int d = 1; d < dim; ++d)
    {
      upper_right[d] = 1;
      repetitions.push_back(1);
    }

  GridGenerator::subdivided_hyper_rectangle(triangulation,
                                            repetitions,
                                            lower_left,
                                            upper_right);
}



int
main()
{
  const int dim = 1;

  Triangulation<dim> triangulation;
  make_two_element_mesh(triangulation);

  hp::FECollection<dim> fe_collection;
  fe_collection.push_back(FE_Nothing<dim>());
  fe_collection.push_back(FE_Q<dim>(1));

  DoFHandler<dim> dof_handler(triangulation);
  dof_handler.begin_active()->set_active_fe_index(1);
  dof_handler.distribute_dofs(fe_collection);

  const auto cell = dof_handler.begin_active();

  std::vector<types::global_dof_index> dof_indices(
    fe_collection[1].dofs_per_cell);

  // Get the dofs at the shared face.
  cell->face(1)->get_dof_indices(dof_indices);
}
--------------------------------------------------------
An error occurred in line <296> of file </home/simon/bin/dealii/include/deal.II/dofs/dof_accessor.templates.h> in function
    static void dealii::internal::DoFAccessorImplementation::Implementation::process_dof_index(const dealii::DoFHandler<dim, spacedim>&, unsigned int, unsigned int, unsigned int, unsigned int, const std::integral_constant<int, structdim>&, GlobalIndexType&, const DoFPProcessor&) [with int dim = 1; int spacedim = 1; int structdim = 0; GlobalIndexType = unsigned int; DoFPProcessor = dealii::internal::DoFAccessorImplementation::Implementation::get_dof_index<1, 1, 0>::<lambda(const auto:8&, auto:9&)>]
The violated condition was: 
    ptr != dof_handler.hp_object_fe_indices[structdim].begin() + dof_handler.hp_object_fe_ptr[structdim][obj_index + 1]
Additional information: 
    You are requesting an active FE index that is not assigned to any of
    the cells connected to this entity.

Stacktrace:
-----------
#0  ./example: void dealii::internal::DoFAccessorImplementation::Implementation::process_dof_index<1, 1, 0, unsigned int, dealii::internal::DoFAccessorImplementation::Implementation::get_dof_index<1, 1, 0>(dealii::DoFHandler<1, 1> const&, unsigned int, unsigned int, unsigned int, unsigned int, std::integral_constant<int, 0> const&)::{lambda(auto:1 const&, auto:2&)#1}>(dealii::DoFHandler<1, 1> const&, unsigned int, unsigned int, unsigned int, unsigned int, std::integral_constant<int, 0> const&, unsigned int&, dealii::internal::DoFAccessorImplementation::Implementation::get_dof_index<1, 1, 0>(dealii::DoFHandler<1, 1> const&, unsigned int, unsigned int, unsigned int, unsigned int, std::integral_constant<int, 0> const&)::{lambda(auto:1 const&, auto:2&)#1} const&)
#1  ./example: unsigned int dealii::internal::DoFAccessorImplementation::Implementation::get_dof_index<1, 1, 0>(dealii::DoFHandler<1, 1> const&, unsigned int, unsigned int, unsigned int, unsigned int, std::integral_constant<int, 0> const&)
#2  ./example: dealii::DoFAccessor<0, 1, 1, false>::get_dof_indices(std::vector<unsigned int, std::allocator<unsigned int> >&, unsigned int) const
#3  ./example: main

@peterrum
Copy link
Member

I think this is rather an issue of the size of dof_indices and that get_dof_indices() does not get the right index. The following code is working:

  {
  std::vector<types::global_dof_index> dof_indices(
    fe_collection[0].dofs_per_face);

    cell->face(1)->get_dof_indices(dof_indices, 0);

    for(const auto i : dof_indices)
      std::cout << i << " ";
    std::cout << std::endl;
  }

  {
  std::vector<types::global_dof_index> dof_indices(
    fe_collection[1].dofs_per_face);

    cell->face(1)->get_dof_indices(dof_indices, 1);

    for(const auto i : dof_indices)
      std::cout << i << " ";
    std::cout << std::endl;
  }

@kronbichler
Copy link
Member

Is there a way to print better assertions that inform the user that the wrong size has been handed in?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants