Segmentation fault in inf_compute_constraints with 2D infinite elements

[AlbertoBittes]

Hello everyone, I have a D-shaped domain and I generated infinite elements on the circular boundary of my domain. I issued a mesh.find_neighbors(); after generating as shown in the examples. I am not doing any AMR. During creation of constraints I am getting a segmentation fault in inf_compute_constraints function, which is called by lagrange_compute_constraints:

template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim, T_radial, T_map>::inf_compute_constraints (DofConstraints & constraints,
                                                           DofMap & dof_map,
                                                           const unsigned int variable_number,
                                                           const Elem * child_elem)
{

  // only constrain elements in 2,3d.
  if (Dim == 1)
    return;

  libmesh_assert(child_elem);

  // only constrain active and ancestor elements
  if (child_elem->subactive())
    return;

  // Before we start to compute anything, lets check if any confinement is needed:
  bool need_constraints=false;
  for (auto child_neighbor : child_elem->neighbor_ptr_range())
    if (child_neighbor->level() < child_elem->level())
      {
        need_constraints = true;
        break;
      }
  if (!need_constraints)
    return;

The neighbor_ptr_range() is returning a nullptr and I am getting a segmentation fault. By compiling libMesh without AMR, the function is skipped and I could complete the constraints creation (BTW, I figured out later that 2D InfFE shape functions are not implemented, is there any workaround or plan to make them available?). Thank you in advance for your help.

[jwpeterson]

Thanks for this bug report, it seems to me that we need to guard this code similarly to the way it is done in elem.h, e.g.

  for (auto n : this->neighbor_ptr_range())
    if (n)
      // do something with n

(BTW, I figured out later that 2D InfFE shape functions are not implemented, is there any workaround or plan to make them available?).

This would be a question for @BalticPinguin, I'm not aware of any ongoing InfFE developments at the moment.

[BalticPinguin]

I am not planning developments with 2D elements..

According to the literature, the situation is more difficult in 2D, because the analytic solution is not of polynomial form, so the numeric solution must be expected to be of worse quality; so I am not sure about quality. Since I never looked at 2D in detail, I am not familiar with the details: If I recall correctly the main change is that you do not have a '1/r'-decay, but instead a '1/sqrt(r)'-decay in the slowest term. In the current implementation, the radial basis is a '1/r'-function, multpilied with a polynomial in 1/r. The simplest (I am not sure if the best) way might be to use the functions as they are, but multiply with 'sqrt(r)'?