Computing the global index values from vertex iterator

[Reza Rastak]

Hi all,

I am maintaining a code containing a function to obtain the global dof numbers of a vertex based on a given vertex iterator (without any use of the cell iterator). The implementation for this function and a sample program using it is shown below:

#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <iostream>
using namespace dealii;
template <int dim>
std::vector<types::global_dof_index> get_vertex_dofs(
  const typename Triangulation<dim>::active_vertex_iterator &vertex,
  const DoFHandler<dim> &                                    dof_handler)
{
  DoFAccessor<0, DoFHandler<dim>, false> vertex_dofs(
      &(dof_handler.get_triangulation()),
      vertex->level(),
      vertex->index(),
      &dof_handler);
  const auto              n_dofs = dof_handler.get_fe().dofs_per_vertex;
  std::vector<types::global_dof_index> dofs(n_dofs);
  for (unsigned int i = 0; i < n_dofs; ++i)
  {
    dofs[i] = vertex_dofs.vertex_dof_index(0, i);
  }
  return dofs;
}

int main() {
  Triangulation<2> triangulation;
  GridGenerator::hyper_cube(triangulation);
  DoFHandler<2> dof_handler(triangulation);
  const FESystem<2> finite_element(FE_Q<2>(1), 2);
  dof_handler.distribute_dofs(finite_element);
  auto vertex = triangulation.begin_active_vertex();
  auto vertex_end = triangulation.end_vertex();
  for (; vertex != vertex_end; ++vertex)
  {
    const auto dofs = get_vertex_dofs(vertex, dof_handler);
    std::cout << dofs[0] << ", " << dofs[1] << std::endl;
  }

  return 0;

}

I am able to run this code with dealii 9.2.0 and see the following output:

0, 1
2, 3
4, 5
6, 7

However, when I try to run it with a newer versions of dealii, the code fails to compile. For example, if I compile the aforementioned code with the latest dealii (commit 855afe3f4c), I get an error because the type DoFAccessor<0, DoFHandler<dim>, false> does not exist anymore. Therefore, I changed the type to DoFAccessor<0, dim, dim, false>. However, the code still fails to compile and shows this error to me:

In file included from /home/reza_rastak/dealii/include/deal.II/dofs/dof_accessor.h:2135,
                 from /home/reza_rastak/dealii/include/deal.II/dofs/dof_handler.h:32,
                 from /home/reza_rastak/dealsims/dealii-examples/vertex-dof/vertex_dof.cc:3:
/home/reza_rastak/dealii/include/deal.II/dofs/dof_accessor.templates.h: In instantiation of 'unsigned int dealii::DoFAccessor<structdim, dim, spacedim, lda>::nth_active_fe_index(unsigned int) const [with int structdim = 0; int dim = 2; int spacedim = 2; bool level_dof_access = false]':
/home/reza_rastak/dealii/include/deal.II/dofs/dof_accessor.templates.h:1544:17:   required from 'dealii::types::global_dof_index dealii::DoFAccessor<structdim, dim, spacedim, lda>::vertex_dof_index(unsigned int, unsigned int, unsigned int) const [with int structdim = 0; int dim = 2; int spacedim = 2; bool level_dof_access = false; dealii::types::global_dof_index = unsigned int]'
/home/reza_rastak/dealsims/dealii-examples/vertex-dof/vertex_dof.cc:27:43:   required from 'std::vector<unsigned int> get_vertex_dofs(const typename dealii::Triangulation<dim, dim>::active_vertex_iterator&, const dealii::DoFHandler<dim>&) [with int dim = 2; typename dealii::Triangulation<dim, dim>::active_vertex_iterator =dealii::TriaActiveIterator<dealii::TriaAccessor<0, 2, 2> >]'
/home/reza_rastak/dealsims/dealii-examples/vertex-dof/vertex_dof.cc:43:58:   required from here
/home/reza_rastak/dealii/include/deal.II/dofs/dof_accessor.templates.h:1488:31: error: 'const class dealii::DoFAccessor<0, 2, 2, false>' has no member named 'present_index'
 1488 |                         this->present_index,
      |                         ~~~~~~^~~~~~~~~~~~~

I dug into the deal.ii codebase, and it seems the error is caused because TriaAccessor<0, dim, spacedim> is not a subclass of TriaAccessorBase<0, dim, spacedim> and thus it does not have a protected member named present_index. However, this protected member is needed within many of the member functions of the DoFAccessor class.

Could you help me to solve this compile issue and keep this function in my codebase?

Thanks,

Reza Rastak

[Peter Munch]

That should fix the issue: https://github.com/dealii/dealii/pull/13429

Peter