BasicRayMeshIntersection.cc

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
/*---------------------------------------------------------------------------*/
/* BasicRayMeshIntersection.cc                                 (C) 2000-2026 */
/*                                                                           */
/* Basic service for calculating intersection between segments and mesh.     */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#include "arcane/utils/Array.h"
#include "arcane/utils/Real3.h"
#include "arcane/utils/StringBuilder.h"
#include "arcane/utils/TraceAccessor.h"
#include "arcane/utils/NotImplementedException.h"
 
#include "arcane/core/ISubDomain.h"
#include "arcane/core/IMesh.h"
#include "arcane/core/IItemFamily.h"
#include "arcane/core/Item.h"
#include "arcane/core/BasicService.h"
#include "arcane/core/FactoryService.h"
#include "arcane/core/ItemPrinter.h"
#include "arcane/core/ItemGroup.h"
#include "arcane/core/VariableTypes.h"
#include "arcane/core/IRayMeshIntersection.h"
#include "arcane/core/IParallelMng.h"
#include "arcane/core/IParticleFamily.h"
 
#include "arcane_packages.h"
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class RayTriangle3DIntersection
: public TraceAccessor
{
 public:
 
  RayTriangle3DIntersection(ITraceMng* tm)
  : TraceAccessor(tm)
  {
  }

  void setRays(Real3ConstArrayView origins, Real3ConstArrayView directions)
  {
    m_rays_origin = origins;
    m_rays_direction = directions;
  }

  void setTriangles(Real3ConstArrayView coordinates, Int32ConstArrayView indexes)
  {
    m_triangles_coordinates = coordinates;
    m_triangles_indexes = indexes;
  }

  void compute();
  RealConstArrayView distances()
  {
    return m_distances;
  }

  Int32ConstArrayView intersectedTriangleIndexes()
  {
    return m_intersected_triangle_indexes;
  }
 
  Real checkIntersection(Real3 origin, Real3 direction, Real3 p0, Real3 p1, Real3 p2);
 
  static bool checkBoundingBox(Real3 origin, Real3 direction, Real3 box_min, Real3 box_max);
 
 private:
 
  void _compute2(Int32 triangle_id, Real3 p0, Real3 p1, Real3 p2);
 
 private:
 
  Real3ConstArrayView m_rays_origin;
  Real3ConstArrayView m_rays_direction;
  Real3ConstArrayView m_triangles_coordinates;
  Int32ConstArrayView m_triangles_indexes;
 
  RealUniqueArray m_distances;
  Int32UniqueArray m_intersected_triangle_indexes;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void RayTriangle3DIntersection::
compute()
{
  Integer nb_ray = m_rays_origin.size();
  m_distances.resize(nb_ray);
  m_distances.fill(1.0e100);
  m_intersected_triangle_indexes.resize(nb_ray);
  m_intersected_triangle_indexes.fill(-1);
  Integer nb_index = m_triangles_indexes.size();
  if ((nb_index % 3) != 0)
    throw FatalErrorException(A_FUNCINFO, "bad triangle_index count (({0}) % 3)!=0");
  Integer nb_triangle = nb_index / 3;
  info() << "COMPUTE RAY INTERSECTION nb_ray=" << nb_ray
         << " nb_triangle=" << nb_triangle;
 
  for (Integer i = 0; i < nb_triangle; ++i) {
    Real3 p0 = m_triangles_coordinates[m_triangles_indexes[(i * 3)]];
    Real3 p1 = m_triangles_coordinates[m_triangles_indexes[(i * 3) + 1]];
    Real3 p2 = m_triangles_coordinates[m_triangles_indexes[(i * 3) + 2]];
    _compute2(i, p0, p1, p2);
  }
 
  for (Integer i = 0; i < nb_ray; ++i) {
    Int32 tid = m_intersected_triangle_indexes[i];
    if (tid == (-1))
      m_distances[i] = -1.0;
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void RayTriangle3DIntersection::
_compute2(Int32 triangle_id, Real3 p0, Real3 p1, Real3 p2)
{
  Integer nb_ray = m_rays_origin.size();
  for (Integer i = 0; i < nb_ray; ++i) {
    Real t = checkIntersection(m_rays_origin[i], m_rays_direction[i], p0, p1, p2);
    if (t >= 0.0) {
      info() << "Segment " << i << " intersect triangle " << triangle_id << " T=" << t;
      if (t < m_distances[i]) {
        m_distances[i] = t;
        m_intersected_triangle_indexes[i] = triangle_id;
        info() << "Get this triangle";
      }
    }
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

Real RayTriangle3DIntersection::
checkIntersection(Real3 origin,
                  Real3 direction,
                  Real3 p0, Real3 p1, Real3 p2)
{
  // This routine is inspired by the code from the WildMagic library
  // whose license is below.
 
  // Wild Magic Source Code
  // David Eberly
  // http://www.geometrictools.com
  // Copyright (c) 1998-2009
  //
  // This library is free software; you can redistribute it and/or modify it
  // under the terms of the GNU Lesser General Public License as published by
  // the Free Software Foundation; either version 2.1 of the License, or (at
  // your option) any later version.  The license is available for reading at
  // either of the locations:
  //     http://www.gnu.org/copyleft/lgpl.html
  //     http://www.geometrictools.com/License/WildMagicLicense.pdf
  //
 
  Real ZERO_TOLERANCE = 1.0e-14;
  // compute the offset origin, edges, and normal
  Real3 kDiff = origin - p0;
  Real3 kEdge1 = p1 - p0;
  Real3 kEdge2 = p2 - p0;
  Real3 kNormal = math::vecMul(kEdge1, kEdge2);
 
  Real fDdN = math::dot(direction, kNormal);
  Real fSign;
  if (fDdN > ZERO_TOLERANCE) {
    fSign = (Real)1.0;
  }
  else if (fDdN < -ZERO_TOLERANCE) {
    fSign = (Real)-1.0;
    fDdN = -fDdN;
  }
  else {
    // Segment and triangle are parallel, call it a "no intersection"
    // even if the segment does intersect.
    return (-1.0);
  }
 
  Real fDdQxE2 = fSign * math::dot(direction, math::vecMul(kDiff, kEdge2));
  if (fDdQxE2 >= (Real)0.0) {
    Real fDdE1xQ = fSign * math::dot(direction, math::vecMul(kEdge1, kDiff));
    if (fDdE1xQ >= (Real)0.0) {
      Real diff = fDdN - (fDdQxE2 + fDdE1xQ);
      // Epsilon is used if the segment passes through the face near
      // one of the face edges. In this case, we consider there
      // to be an intersection. Without this, due to numerical rounding,
      // a segment could pass through the mesh by going through an
      // edge between two faces without intersecting them.
      if (diff >= 0.0 || diff > -ZERO_TOLERANCE) {
        // line intersects triangle, check if segment does
        Real fQdN = -fSign * math::dot(kDiff, kNormal);
 
        // Since this is a half-line and not a segment,
        // there is always an intersection if \a t is positive.
#if 0
        Real fExtDdN = m_pkSegment->Extent*fDdN;
        if (-fExtDdN <= fQdN && fQdN <= fExtDdN)
        {
          // segment intersects triangle
          return true;
        }
#endif
        Real fInv = ((Real)1.0) / fDdN;
        Real t = fQdN * fInv;
        //Real b1 = fDdQxE2*fInv;
        //Real b2 = fDdE1xQ*fInv;
        //Real b0 = (Real)1.0 - b1 - b2;
        //   P = origin + t*direction = b0*V0 + b1*V1 + b2*V2
        //info() << "** INTERSECT POS=" << b0 << " y=" << b1 << " z=" << b2 << " T=" << t;
        //Real intersect_p = origin + t*direction;
        return t;
 
        // else: |t| > extent, no intersection
      }
      // else: b1+b2 > 1, no intersection
    }
    // else: b2 < 0, no intersection
  }
  // else: b1 < 0, no intersection
 
  return (-1.0);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
bool RayTriangle3DIntersection::
checkBoundingBox(Real3 origin, Real3 direction, Real3 box_min, Real3 box_max)
{
  Real3 box_center = (box_max + box_min) * 0.5;
  Real afWdU[3], afAWdU[3], afDdU[3], afADdU[3], afAWxDdU[3], fRhs;
 
  Real3 kDiff = origin - box_center;
  Real3 Axis[3];
  //Axis[0].x = box_max.x - box_min.x;
  //Axis[1].y = box_max.y - box_min.y;
  //Axis[2].z = box_max.z - box_min.z;
  Axis[0] = Real3(1.0, 0.0, 0.0);
  Axis[1] = Real3(0.0, 1.0, 0.0);
  Axis[2] = Real3(0.0, 0.0, 1.0);
  Real Extent[3];
  Extent[0] = box_max.x - box_min.x;
  Extent[1] = box_max.y - box_min.y;
  Extent[2] = box_max.z - box_min.z;
 
  afWdU[0] = math::dot(direction, Axis[0]);
  afAWdU[0] = math::abs(afWdU[0]);
  afDdU[0] = math::dot(kDiff, Axis[0]);
  afADdU[0] = math::abs(afDdU[0]);
  if (afADdU[0] > Extent[0] && afDdU[0] * afWdU[0] >= (Real)0.0) {
    return false;
  }
 
  afWdU[1] = math::dot(direction, Axis[1]);
  afAWdU[1] = math::abs(afWdU[1]);
  afDdU[1] = math::dot(kDiff, Axis[1]);
  afADdU[1] = math::abs(afDdU[1]);
  if (afADdU[1] > Extent[1] && afDdU[1] * afWdU[1] >= (Real)0.0) {
    return false;
  }
 
  afWdU[2] = math::dot(direction, Axis[2]);
  afAWdU[2] = math::abs(afWdU[2]);
  afDdU[2] = math::dot(kDiff, Axis[2]);
  afADdU[2] = math::abs(afDdU[2]);
  if (afADdU[2] > Extent[2] && afDdU[2] * afWdU[2] >= (Real)0.0) {
    return false;
  }
 
  Real3 kWxD = math::vecMul(direction, kDiff);
 
  afAWxDdU[0] = math::abs(math::dot(kWxD, Axis[0]));
  fRhs = Extent[1] * afAWdU[2] + Extent[2] * afAWdU[1];
  if (afAWxDdU[0] > fRhs) {
    return false;
  }
 
  afAWxDdU[1] = math::abs(math::dot(kWxD, Axis[1]));
  fRhs = Extent[0] * afAWdU[2] + Extent[2] * afAWdU[0];
  if (afAWxDdU[1] > fRhs) {
    return false;
  }
 
  afAWxDdU[2] = math::abs(math::dot(kWxD, Axis[2]));
  fRhs = Extent[0] * afAWdU[1] + Extent[1] * afAWdU[0];
  if (afAWxDdU[2] > fRhs) {
    return false;
  }
 
  return true;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
class BasicRayFaceIntersector
: public IRayFaceIntersector
{
 public:
 
  BasicRayFaceIntersector(ITraceMng* tm)
  : m_triangle_intersector(tm)
  {
  }
 
 public:
 
  bool computeIntersection(Real3 origin, Real3 direction,
                           Int32 orig_face_local_id,
                           Int32 face_local_id,
                           Real3ConstArrayView face_nodes,
                           Int32* user_value,
                           Real* distance, Real3* position)
  {
    ARCANE_UNUSED(orig_face_local_id);
    ARCANE_UNUSED(face_local_id);
 
    Integer nb_node = face_nodes.size();
    switch (nb_node) {
    case 4: {
      Real3 center = (face_nodes[0] + face_nodes[1] + face_nodes[2] + face_nodes[3]) / 4.0;
      Real d0 = m_triangle_intersector.checkIntersection(origin, direction, center, face_nodes[0], face_nodes[1]);
      Real d1 = m_triangle_intersector.checkIntersection(origin, direction, center, face_nodes[1], face_nodes[2]);
      Real d2 = m_triangle_intersector.checkIntersection(origin, direction, center, face_nodes[2], face_nodes[3]);
      Real d3 = m_triangle_intersector.checkIntersection(origin, direction, center, face_nodes[3], face_nodes[0]);
      Real d = 1.0e100;
      bool found = false;
      if (d0 >= 0.0 && d0 < d) {
        found = true;
        d = d0;
      }
      if (d1 >= 0.0 && d1 < d) {
        found = true;
        d = d1;
      }
      if (d2 >= 0.0 && d2 < d) {
        found = true;
        d = d2;
      }
      if (d3 >= 0.0 && d3 < d) {
        found = true;
        d = d3;
      }
      if (!found)
        d = -1.0;
      *distance = d;
      *position = origin + d * direction;
      *user_value = 1;
      return found;
    } break;
    default:
      throw NotImplementedException(A_FUNCINFO, "only quad face is implemented");
    }
  }
 
 public:
 
  RayTriangle3DIntersection m_triangle_intersector;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class BasicRayMeshIntersection
: public BasicService
, public IRayMeshIntersection
{
 public:
 
  BasicRayMeshIntersection(const ServiceBuildInfo& sbi);
  virtual ~BasicRayMeshIntersection() {}
 
 public:
 
  virtual void build() {}
  virtual void compute(Real3ConstArrayView segments_position,
                       Real3ConstArrayView segments_direction,
                       Int32ConstArrayView segments_orig_face,
                       Int32ArrayView user_values,
                       Real3ArrayView intersections,
                       RealArrayView distances,
                       Int32ArrayView faces_local_id);
 
  virtual void compute(IItemFamily* ray_family,
                       VariableParticleReal3& rays_position,
                       VariableParticleReal3& rays_direction,
                       VariableParticleInt32& rays_orig_face,
                       VariableParticleInt32& user_values,
                       VariableParticleReal3& intersections,
                       VariableParticleReal& distances,
                       VariableParticleInt32& rays_face);
 
  virtual void setFaceIntersector(IRayFaceIntersector* intersector)
  {
    m_face_intersector = intersector;
  }
  virtual IRayFaceIntersector* faceIntersector()
  {
    return m_face_intersector;
  }
 
 public:
 
  inline void _checkBoundingBox(Real3 p, Real3* ARCANE_RESTRICT min_bounding_box,
                                Real3* ARCANE_RESTRICT max_bounding_box)
  {
    if (p.x < min_bounding_box->x)
      min_bounding_box->x = p.x;
    if (p.y < min_bounding_box->y)
      min_bounding_box->y = p.y;
    if (p.z < min_bounding_box->z)
      min_bounding_box->z = p.z;
 
    if (p.x > max_bounding_box->x)
      max_bounding_box->x = p.x;
    if (p.y > max_bounding_box->y)
      max_bounding_box->y = p.y;
    if (p.z > max_bounding_box->z)
      max_bounding_box->z = p.z;
  }
 
 private:
 
  IRayFaceIntersector* m_face_intersector;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
BasicRayMeshIntersection::
BasicRayMeshIntersection(const ServiceBuildInfo& sbi)
: BasicService(sbi)
, m_face_intersector(0)
{
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void BasicRayMeshIntersection::
compute(Real3ConstArrayView segments_position,
        Real3ConstArrayView segments_direction,
        Int32ConstArrayView segments_orig_face,
        Int32ArrayView user_values,
        Real3ArrayView segments_intersection,
        RealArrayView segments_distance,
        Int32ArrayView faces_local_id)
{
  IMesh* mesh = this->mesh();
 
  bool is_3d = mesh->dimension() == 3;
  info() << "COMPUTE INTERSECTION!!";
  FaceGroup outer_faces = mesh->outerFaces();
  Integer nb_face = outer_faces.size();
  Integer nb_segment = segments_position.size();
  info() << "NB OUTER FACE=" << nb_face << " NB_SEGMENT=" << nb_segment;
  VariableNodeReal3 nodes_coordinates(mesh->nodesCoordinates());
  const Real max_value = 1.0e100;
  const Real3 max_bb(max_value, max_value, max_value);
  const Real3 min_bb(-max_value, -max_value, -max_value);
 
  Real3 mesh_min_bounding_box(max_bb);
  Real3 mesh_max_bounding_box(min_bb);
  Integer max_face_local_id = mesh->faceFamily()->maxLocalId();
  Real3UniqueArray faces_min_bb(max_face_local_id);
  Real3UniqueArray faces_max_bb(max_face_local_id);
 
  // Calculate the bounding boxes
  {
    ENUMERATE_FACE (iface, outer_faces) {
      Int32 lid = iface.localId();
      Real3 face_max_bb(max_bb);
      Real3 face_min_bb(min_bb);
      for (NodeEnumerator inode((*iface).nodes()); inode.hasNext(); ++inode)
        _checkBoundingBox(nodes_coordinates[inode], &face_min_bb, &face_max_bb);
      _checkBoundingBox(face_min_bb, &mesh_min_bounding_box, &mesh_max_bounding_box);
      _checkBoundingBox(face_max_bb, &mesh_min_bounding_box, &mesh_max_bounding_box);
      //TODO: maybe add an epsilon around the BB to avoid
      // rounding errors in determining the intersection
      faces_min_bb[lid] = face_min_bb;
      faces_max_bb[lid] = face_max_bb;
    }
  }
 
  //RealArray segments_distance;
  segments_distance.fill(max_value);
  faces_local_id.fill(NULL_ITEM_LOCAL_ID);
 
  if (!m_face_intersector)
    m_face_intersector = new BasicRayFaceIntersector(traceMng());
 
  Real3UniqueArray face_nodes;
  for (Integer i = 0; i < nb_segment; ++i) {
    Real3 position = segments_position[i];
    Real3 direction = segments_direction[i];
    Int32 orig_face_local_id = segments_orig_face[i];
    Real3 intersection;
    Real distance = 1.0e100;
    Int32 min_face_local_id = NULL_ITEM_LOCAL_ID;
    Int32 user_value = 0;
    ENUMERATE_FACE (iface, outer_faces) {
      const Face& face = *iface;
      // Only process its own faces
      if (!face.isOwn())
        continue;
      Int32 lid = face.localId();
      // In 3D, check for intersection with the bounding box.
      // If there is none, there is no need to go further.
      bool is_bb = true;
      if (is_3d)
        is_bb = RayTriangle3DIntersection::checkBoundingBox(position, direction, faces_min_bb[lid], faces_max_bb[lid]);
      if (!is_bb)
        continue;
      Integer nb_node = face.nbNode();
      face_nodes.resize(nb_node);
      for (NodeEnumerator inode(face.nodes()); inode.hasNext(); ++inode)
        face_nodes[inode.index()] = nodes_coordinates[inode];
      Real d = 0.0;
      Real3 local_intersection;
      Int32 uv = 0;
      bool is_found = m_face_intersector->computeIntersection(position, direction, orig_face_local_id, lid,
                                                              face_nodes, &uv, &d, &local_intersection);
      //if (i==15){
      //  for( Integer z=0; z<nb_node; ++z )
      //    info() << "FACE NODE lid=" << lid << " I=" << i << " v=" << face_nodes[z] << " d=" << d;
      //}
      //if (!is_bb)
      //info() << "CHECK INTERSECTION: is_bb=" << is_bb << " " << is_found << '\n';
      if (is_found && !is_bb)
        fatal() << "Intersection found but no bounding box intersection";
      if (is_found && d < distance) {
        distance = d;
        min_face_local_id = lid;
        intersection = local_intersection;
        user_value = uv;
      }
    }
    if (min_face_local_id == NULL_ITEM_LOCAL_ID)
      distance = -1.0;
    segments_distance[i] = distance;
    segments_intersection[i] = intersection;
    faces_local_id[i] = min_face_local_id;
    user_values[i] = user_value;
  }
}
 
struct FoundInfo
{
  Int32 contrib_owner;
  Int32 owner;
  Int32 local_id;
  Int32 face_local_id;
  Int32 orig_face_local_id;
  Int32 user_value;
  Real distance;
  Real3POD intersection;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void BasicRayMeshIntersection::
compute(IItemFamily* ray_family,
        VariableParticleReal3& rays_position,
        VariableParticleReal3& rays_direction,
        VariableParticleInt32& rays_orig_face,
        VariableParticleInt32& rays_user_value,
        VariableParticleReal3& intersections,
        VariableParticleReal& distances,
        VariableParticleInt32& rays_face)
{
  // NOTE: rays_orig_face is not used.
 
  IMesh* mesh = ray_family->mesh();
  IParallelMng* pm = mesh->parallelMng();
  Int32 my_rank = pm->commRank();
  // Assume that the rays are compacted
  Integer nb_local_ray = ray_family->allItems().size();
  //if (nb_local_ray!=ray_family->maxLocalId()){
  //fatal() << "The ray family must be compacted nb=" << nb_local_ray
  //          << " max_id=" << ray_family->maxLocalId();
  //}
  Integer global_nb_ray = pm->reduce(Parallel::ReduceSum, nb_local_ray);
  if (global_nb_ray == 0)
    return;
  info() << "LOCAL_NB_RAY=" << nb_local_ray << " GLOBAL=" << global_nb_ray;
 
  // Copy input information into a local array.
  // This is always necessary because particles are not necessarily
  // compacted and there may be holes in the numbering.
  Real3UniqueArray local_positions(nb_local_ray);
  Real3UniqueArray local_directions(nb_local_ray);
 
  ENUMERATE_PARTICLE (iitem, ray_family->allItems()) {
    Integer index = iitem.index();
    //local_ids[index] = iitem.localId();
    //unique_ids[index] = (*iitem).uniqueId();
    local_positions[index] = rays_position[iitem];
    local_directions[index] = rays_direction[iitem];
  }
 
  if (pm->isParallel()) {
 
    Int32UniqueArray local_ids(nb_local_ray);
    Int64UniqueArray unique_ids(nb_local_ray);
    ENUMERATE_ITEM (iitem, ray_family->allItems()) {
      Integer index = iitem.index();
      local_ids[index] = iitem.localId();
      unique_ids[index] = (*iitem).uniqueId();
    }
 
    // In parallel, currently retrieve rays from other processors
 
    Real3UniqueArray all_positions;
    pm->allGatherVariable(local_positions, all_positions);
    Real3UniqueArray all_directions;
    pm->allGatherVariable(local_directions, all_directions);
 
    Int32UniqueArray local_owners(nb_local_ray);
    local_owners.fill(my_rank);
 
    Int32UniqueArray all_owners;
    pm->allGatherVariable(local_owners, all_owners);
    Int32UniqueArray all_local_ids;
    pm->allGatherVariable(local_ids, all_local_ids);
    Int64UniqueArray all_unique_ids;
    pm->allGatherVariable(unique_ids, all_unique_ids);
 
    Int32UniqueArray all_user_values(global_nb_ray);
    Int32UniqueArray all_orig_faces(global_nb_ray, NULL_ITEM_LOCAL_ID);
    RealUniqueArray all_distances(global_nb_ray);
    Real3UniqueArray all_intersections(global_nb_ray);
    Int32UniqueArray all_faces(global_nb_ray);
    // Set faces that I am not the owner of to (-1)
    for (Integer z = 0, zn = all_orig_faces.size(); z < zn; ++z) {
      if (all_owners[z] != my_rank)
        all_orig_faces[z] = NULL_ITEM_LOCAL_ID;
    }
    compute(all_positions, all_directions, all_orig_faces, all_user_values, all_intersections, all_distances, all_faces);
    /*for( Integer i=0; i<global_nb_ray; ++i ){
      info() << "RAY I=" << i << " uid=" << all_unique_ids[i] << " lid=" << all_local_ids[i]
             << " owner=" << all_owners[i] << " position=" << all_positions[i]
             << " direction=" << all_directions[i] << " distance=" << all_distances[i]
             << " face=" << all_faces[i];
             }*/
    {
      UniqueArray<FoundInfo> founds_info;
 
      for (Integer i = 0; i < global_nb_ray; ++i) {
        Int32 face_lid = all_faces[i];
        if (face_lid == NULL_ITEM_LOCAL_ID)
          continue;
        FoundInfo fi;
        fi.contrib_owner = my_rank;
        fi.owner = all_owners[i];
        fi.local_id = all_local_ids[i];
        fi.face_local_id = all_faces[i];
        fi.orig_face_local_id = all_orig_faces[i];
        fi.user_value = all_user_values[i];
        fi.intersection = all_intersections[i];
        fi.distance = all_distances[i];
        founds_info.add(fi);
      }
      info() << "NB_FOUND=" << founds_info.size();
      //Array<FoundInfo> global_founds_info;
      ByteUniqueArray global_founds_info_bytes;
      Integer finfo_byte_size = arcaneCheckArraySize(founds_info.size() * sizeof(FoundInfo));
      ByteConstArrayView local_fi(finfo_byte_size, (const Byte*)founds_info.data());
      pm->allGatherVariable(local_fi, global_founds_info_bytes);
      Integer gfi_byte_size = arcaneCheckArraySize(global_founds_info_bytes.size() / sizeof(FoundInfo));
      ConstArrayView<FoundInfo> global_founds_info(gfi_byte_size, (const FoundInfo*)global_founds_info_bytes.data());
      Integer nb_total_found = global_founds_info.size();
      ParticleInfoListView rays(ray_family);
      rays_face.fill(NULL_ITEM_LOCAL_ID);
      distances.fill(0.0);
      VariableItemInt32& rays_new_owner = ray_family->itemsNewOwner();
      rays_new_owner.fill(my_rank);
      for (Integer i = 0; i < nb_total_found; ++i) {
        const FoundInfo& fi = global_founds_info[i];
        Int32 owner = fi.owner;
        // Only process its own rays
        if (owner != my_rank)
          continue;
        Int32 local_id = fi.local_id;
        Particle ray = rays[local_id];
        Real distance = fi.distance;
        if ((rays_face[ray] == NULL_ITEM_LOCAL_ID) || distance < distances[ray]) {
          distances[ray] = distance;
          intersections[ray] = Real3(fi.intersection.x, fi.intersection.y, fi.intersection.z);
          rays_face[ray] = fi.face_local_id;
          rays_user_value[ray] = fi.user_value;
          rays_new_owner[ray] = fi.contrib_owner;
          //info() << "SET RAY uid=" << ray.uniqueId() << " distance=" << distance << " new_owner=" << fi.contrib_owner;
        }
        /*info() << "GLOBAL RAY I=" << i << " uid=" << ray.uniqueId() << " lid=" << fi.local_id
               << " owner=" << fi.owner
               << " distance=" << fi.distance
               << " contrib_owner=" << fi.contrib_owner
               << " face=" << fi.face_local_id
               << " new_owner=" << rays_new_owner[ray];*/
      }
      ray_family->toParticleFamily()->exchangeParticles();
    }
  }
  else {
 
    Int32UniqueArray orig_faces(nb_local_ray);
    Int32UniqueArray user_values(nb_local_ray);
 
    ENUMERATE_PARTICLE (iitem, ray_family->allItems()) {
      Integer index = iitem.index();
      orig_faces[index] = rays_orig_face[iitem];
      user_values[index] = rays_user_value[iitem];
    }
 
    Real3UniqueArray out_intersections(nb_local_ray);
    RealUniqueArray out_distances(nb_local_ray);
    Int32UniqueArray out_faces(nb_local_ray);
 
    compute(local_positions, local_directions, orig_faces, user_values, out_intersections, out_distances, out_faces);
 
    // Copy values back to the corresponding variables in the output.
    ENUMERATE_PARTICLE (iitem, ray_family->allItems()) {
      Integer index = iitem.index();
      intersections[iitem] = out_intersections[index];
      distances[iitem] = out_distances[index];
      rays_face[iitem] = out_faces[index];
    }
  }
 
  // For testing, write the rays and their impact point.
  {
    Integer nb_new_ray = ray_family->nbItem();
    RealUniqueArray local_distances(nb_new_ray);
    local_directions.resize(nb_new_ray);
    local_positions.resize(nb_new_ray);
    ENUMERATE_PARTICLE (iitem, ray_family->allItems()) {
      Integer index = iitem.index();
      local_distances[index] = distances[iitem];
      local_directions[index] = rays_direction[iitem];
      local_positions[index] = rays_position[iitem];
    }
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
ARCANE_REGISTER_SUB_DOMAIN_FACTORY(BasicRayMeshIntersection, IRayMeshIntersection,
                                   BasicRayMeshIntersection);
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // End namespace Arcane
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/