Avoid dying on failed contact point searches

[GiudGiud]

Bug Description

The Newton search in FindContactPoint can try to call reinit on an element side with a point that is way outside of the element.
This causes a libmesh error and crashes the initialization of the simulation

A warning instead of an error could be good enough to get the users to work more on their mesh

Steps to Reproduce

Run the input from this discussion (not_working.i)
#28091 (reply in thread)

or make a similar system
Input

[Mesh]
  [base]
    type = FileMeshGenerator
    file = three_cylinders.e
  []
  [diag]
    type = MeshDiagnosticsGenerator
    input = base
    examine_element_overlap = WARNING
    examine_element_types = WARNING
    examine_element_volumes = WARNING
    examine_non_conformality = WARNING
    examine_nonplanar_sides = INFO
    examine_sidesets_orientation = WARNING
    search_for_adaptivity_nonconformality = WARNING
    check_local_jacobian = WARNING
  []
  patch_update_strategy = iteration
  patch_size = 1000
[]

[GlobalParams]
  volumetric_locking_correction = true
  displacements = 'disp_x disp_y disp_z'
[]

[Modules/TensorMechanics/Master]
  [all]
    add_variables = true
    strain = FINITE
    block = 'cylinder_top cylinder_center cylinder_lower'
    use_automatic_differentiation = false
    use_finite_deform_jacobian = true
  []
[]

[BCs]
  [top_z_d]
    type = FunctionDirichletBC
    variable = disp_z
    boundary = 'cylinder_upper_top'
    function = '-t'
  []
  [top_y_d]
    type = DirichletBC
    variable = disp_y
    boundary = 'cylinder_upper_top'
    value = 0.0
  []
  [top_x_d]
    type = DirichletBC
    variable = disp_x
    boundary = 'cylinder_upper_top'
    value = 0.0
  []
  [bottom_y_d]
    type = DirichletBC
    variable = disp_y
    boundary = 'cylinder_lower_bottom'
    value = 0.0
  []
  [bottom_x_d]
    type = DirichletBC
    variable = disp_x
    boundary = 'cylinder_lower_bottom'
    value = 0.0
  []
  [bottom_z_d]
    type = DirichletBC
    variable = disp_z
    boundary = 'cylinder_lower_bottom'
    value = 0.0
  []
[]

[Materials]
  [center]
    type = ComputeIsotropicElasticityTensor
    block = 'cylinder_center'
    youngs_modulus = 10.0e3
    poissons_ratio = 0.2
  []
  [top]
    type = ComputeIsotropicElasticityTensor
    block = 'cylinder_top cylinder_lower'
    youngs_modulus = 210.0e3
    poissons_ratio = 0.2
  []
  [stress]
    type = ComputeFiniteStrainElasticStress
    block = 'cylinder_top cylinder_center cylinder_lower'
  []
[]

[Preconditioning]
  [SMP]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = 'PJFNK'

  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package -pc_factor_shift_type'
  petsc_options_value = 'lu     superlu_dist NONZERO'

  line_search = 'none'

  l_max_its = 100
  nl_max_its = 120
  end_time = 4.0
  nl_rel_tol = 1.e-6
  nl_abs_tol = 1.e-7
  dtmax = 0.1
  dtmin = 1.e-6
  l_tol = 1e-3
[]

[Outputs]
  file_base = three_cylinders_out
  print_linear_residuals = false
  perf_graph = true
  [exodus]
    type = Exodus
    elemental_as_nodal = true
  []
  [console]
    type = Console
    max_rows = 5
  []
[]

[Contact]
  [top_center_cont]
    primary = 'ss_cylinder_top_lower'
    secondary = 'ss_cylinder_center_upper'
    model = glued
  []
  [lower_center_cont]
    secondary = 'ss_cylinder_center_lower'
    primary = 'ss_cylinder_lower_top'
    model = glued
  []
[]

Mesh is 3 cylinders with completely arbitrary hex8 meshes

Impact

User support needs

[Optional] Diagnostics

No response