class DegenerateMapping : public std::runtime_error
{
public:
  DegenerateMapping(std::string msg="") :
    std::runtime_error( "Problem with mapping or its Jacobian.\n" + msg ) {}
};

    libmesh_assert_equal_to(phi_map.size(), elem_nodes.size());

  auto check_for_degenerate_map =
    [this, elem, p]
    (Real det_J) {
    if (det_J <= jacobian_tolerance)
      {
        // Don't call get_info() recursively if we're already
        // failing.  get_info() calls Elem::volume() which may
        // call FE::reinit() and trigger the same failure again.
        static bool failing = false;
        if (!failing)
          {
            failing = true;
            std::string elem_info = elem->get_info();
            failing = false;
            if (calculate_xyz)
              {
                libmesh_degenerate_mapping_msg
                  ("Jacobian " << jac[p] << " under tolerance " <<
                   jacobian_tolerance << " at point " << xyz[p] <<
                   " in element: " << elem_info);

              {
                // In this case xyz[p] is not defined, so don't
                // try to print it out.
                libmesh_degenerate_mapping_msg
                  ("Jacobian " << jac[p] << " under tolerance " <<
                   jacobian_tolerance << " at point index " << p <<
                   " in element: " << elem_info);

            // We were already failing when we called this, so just
            // stop the current computation and return with
            // incomplete results.
            return;
          }
      }
  };

  switch (dim)
    {

          {
            jac[p] = dxyzdxi_map[p].norm();

            check_for_degenerate_map(jac[p]);

            // The inverse Jacobian entries also come from the
            // generalized inverse of T (see also the 2D element

            const Real det = (g11*g22 - g12*g21);

            check_for_degenerate_map(det);

            const Real inv_det = 1./det;
            jac[p] = std::sqrt(det);

                      dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta)  +
                      dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));

            check_for_degenerate_map(jac[p]);

            JxW[p] = jac[p]*qw[p];

            //  G = [J]^T [J]
            const Real G = dxi*dxi;

            if (secure && G <= 0)
              libmesh_degenerate_mapping_msg
                ("inverse_map found a singular Jacobian " <<
                 " at master point " << p << " in element " <<
                 elem->id());

            const Real det = (G11*G22 - G12*G21);

            if (secure && det == 0)
              libmesh_degenerate_mapping_msg
                ("inverse_map found a singular Jacobian " <<
                 " at master point " << p << " in element " <<
                 elem->id());

                // case we can just return a far away point.
                if (secure)
                  {
                    libmesh_degenerate_mapping_msg(
                      "inverse_map found a singular Jacobian" <<
                      " at master point " << p << " in element " <<
                      elem->id());

                  // "Jacobian becomes singular or negative" sense,
                  // it's at least degenerate in the "straightforward
                  // Newton is failing to invert it" sense.
                  libmesh_degenerate_mapping_msg(
                    "inverse_map Newton FAILED to converge in " <<
                    cnt << " iterations in element " << elem->id() <<
                    " for physical point = " << physical_point <<

  // this captures the case that the point is not (almost) in the direction of the element.:
  // first, project the problem onto the unit sphere:
  Point p_o(p - my_origin);
  pt0_o /= pt0_o.norm();
  pt1_o /= pt1_o.norm();
  pt2_o /= pt2_o.norm();
  p_o /= p_o.norm();

  // now, check if it is in the projected face; by comparing the distance of
  // any point in the element to \p p with the largest distance between this point
  // to any other point in the element.
  if ((p_o - pt0_o).norm_sq() > std::max((pt0_o - pt1_o).norm_sq(), (pt0_o - pt2_o).norm_sq()) ||
      (p_o - pt1_o).norm_sq() > std::max((pt1_o - pt2_o).norm_sq(), (pt1_o - pt0_o).norm_sq()) ||
      (p_o - pt2_o).norm_sq() > std::max((pt2_o - pt0_o).norm_sq(), (pt2_o - pt1_o).norm_sq()) )
    {
      // the physical point is definitely not contained in the element
      return false;
    }

                              if (neighbor->p_level() < my_p_level &&
                                  neighbor->p_refinement_flag() != Elem::REFINE)
                                {
                                  neighbor->set_p_refinement_flag(Elem::REFINE);
                                  level_one_satisfied = false;
                                  compatible_with_coarsening = false;
                                }
                              if (neighbor->p_level() == my_p_level &&

    {
      ray_target = inside - Point(1);
      intersection_distances =
        this->find_ray_intersections(inside, ray_target);
    }

  // I'd make this an assert, but I'm not 100% confident we can't

          libMesh::out << "Maximum EIM error is " << greedy_error << std::endl << std::endl;
        }
      libmesh_catch (const std::exception & e)
        {
          // If we hit an exception when performing the enrichment for the error indicator, then
          // we just continue and skip the error indicator. Otherwise we rethrow the exception.

              update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
            }
          libmesh_catch (const std::exception & e)
            {
              // If we hit an exception when performing the enrichment for the error indicator, then
              // we just continue and skip the error indicator. Otherwise we rethrow the exception.

              update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
            }
          libmesh_catch (const std::exception & e)
            {
              // If we hit an exception when performing the enrichment for the error indicator, then
              // we just continue and skip the error indicator. Otherwise we rethrow the exception.

              update_eim_matrices(/*set_error_indicator*/ exit_on_next_iteration);
            }
          libmesh_catch (const std::exception & e)
            {
              // If we hit an exception when performing the enrichment for the error indicator, then
              // we just continue and skip the error indicator. Otherwise we rethrow the exception.