GeometryKernelSurfaceToolsService.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
//-----------------------------------------------------------------------------
 
#ifdef WIN32
#include <iso646.h>
#endif
 
#include "arcane/corefinement/surfaceutils/geometrykernelsurfacetools/GeometryKernelSurfaceToolsService.h"
/* Author : havep at Wed Aug 27 14:57:32 2008
 * Generated by createNew
 */
 
#include "arcane/corefinement/surfaceutils/geometrykernelsurfacetools/GeometryKernelSurfaceImpl.h"
 
#include "arcane/utils/NotImplementedException.h"
#include "arcane/core/IParallelMng.h"
#include "arcane/core/IMesh.h"
 
#include <GeometryKernel/algorithms/surface/surface-corefinement.h>
#include <GeometryKernel/tools/surface/triangulation-topo-tools.h>
#include <GeometryKernel/datamodel/geometry/vector.h>
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
using namespace Arcane;
 
#include "arcane/corefinement/surfaceutils/geometrykernelsurfacetools/GeometryKernelSurfaceInternalUtils.h"
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void GeometryKernelSurfaceToolsService::
init()
{
  ;
}
 
/*---------------------------------------------------------------------------*/
 
ISurface*
GeometryKernelSurfaceToolsService::
createSurface()
{
  return new GeometryKernelSurfaceImpl();
}
 
/*---------------------------------------------------------------------------*/
 
void GeometryKernelSurfaceToolsService::
setFaceToSurface(ISurface* surface,
                 FaceGroup face_group)
{
  buildFaceGroupSurface(face_group, castSurface(surface));
}
 
/*---------------------------------------------------------------------------*/
 
void GeometryKernelSurfaceToolsService::
computeSurfaceContact(ISurface* s1,
                      ISurface* s2,
                      FaceFaceContactList& contact)
{
  GeometryKernelSurfaceImpl* surface1 = castSurface(s1);
  GeometryKernelSurfaceImpl* surface2 = castSurface(s2);
 
  GK::TriangulationDataStructurePtr gkSurface1 = surface1->m_triangulation;
  GK::TriangulationDataStructurePtr gkSurface2 = surface2->m_triangulation;
 
  // hack to skip unwanted boost pointer issues
  GK::TriangulationDataStructurePtr surfaceA = gkSurface1;
  GK::TriangulationDataStructurePtr surfaceB = gkSurface2;
 
  //   { // GK debug purpose
  //     const bool use_surface_count = false;
  //     const bool force_cpu_id = true;
 
  //     char cpu_suffix[16] = "";
  //     if (force_cpu_id or subDomain()->nbSubDomain() > 1)
  //       sprintf(cpu_suffix,"_cpu=%d",subDomain()->subDomainId());
 
  //     char count_suffix[32] = "";
  //     if (use_surface_count)
  //       {
  //         static Integer surface_count = 0;
  //         sprintf(count_suffix,"_%d",surface_count);
  //         ++surface_count;
  //       }
 
  //     char filename[256];
  //     sprintf(filename,"surfaceA%s%s",cpu_suffix,count_suffix);
  //     ::saveSurface(filename,*surfaceA);
  //     sprintf(filename,"surfaceB%s%s",cpu_suffix,count_suffix);
  //     ::saveSurface(filename,*surfaceB);
  //   }
 
  UniqueArray<Real3> normals(2);
  normals[0] = surface1->m_mean_normal;
  normals[1] = surface2->m_mean_normal;
  //subDomain()->parallelMng()->reduce(Parallel::ReduceSum,normals);
 
#ifndef NO_USER_WARNING
#warning "Work around strange UniqueArray<Real3> reduction"
#endif /* NO_USER_WARNING */
  UniqueArray<Real> __normals(6);
  __normals[0] = normals[0].x;
  __normals[1] = normals[0].y;
  __normals[2] = normals[0].z;
  __normals[3] = normals[1].x;
  __normals[4] = normals[1].y;
  __normals[5] = normals[1].z;
  subDomain()->parallelMng()->reduce(Parallel::ReduceSum, __normals);
  normals[0].x = __normals[0];
  normals[0].y = __normals[1];
  normals[0].z = __normals[2];
  normals[1].x = __normals[3];
  normals[1].y = __normals[4];
  normals[1].z = __normals[5];
 
  // Clear container
  contact.clear();
 
  // No premature return before this point since reduction needs all processors
  if (surfaceA->numberOfFaces() == 0) // non empty master surface
    return;
 
  typedef GK::SurfaceCorefinement::IntersectionIterator Iter;
  typedef GK::SurfaceCorefinement::IntersectionHandle Handle;
  long requested_info = GK::SurfaceCorefinement::NORMAL | GK::SurfaceCorefinement::CENTRE;
#ifndef NO_USER_WARNING
#warning "global plane computation disabled"
#endif /* NO_USER_WARNING */
  GK::SurfaceCorefinement s(surfaceA, surfaceB,
                            convertGKVector(normals[0]), convertGKVector(normals[1]),
                            requested_info, options()->areaEpsilon());
 
  const Integer voidIndex = s.indexVoid();
 
  const Array<Face>& face_arrayA = surface1->m_face_array;
  const Array<Face>& face_arrayB = surface2->m_face_array;
  const Array<bool>& face_reorientA = surface1->m_face_reorient;
  const Array<bool>& face_reorientB = surface2->m_face_reorient;
 
  Integer count = 0;
 
  // voidIndex on indexA is handled separately
  for (Integer indexA = s.beginA(); indexA < s.endA(); /*incr inside*/) {
    ARCANE_ASSERT((indexA != voidIndex), ("indexA cannot be voidIndex by documentation"));
    Face faceA = face_arrayA[indexA];
    Integer split_sizeA = faceA.nbNode() - 2; // assume what split method has been used
 
    // Table of current positions of sub-faces B associated with this face A
    UniqueArray<Iter> jB(split_sizeA);
    UniqueArray<Iter> endB(split_sizeA);
    bool canContinue = false;
    for (Integer iA = 0; iA < split_sizeA; ++iA) {
      jB[iA] = s.beginB(indexA + iA);
      endB[iA] = s.endB(indexA + iA);
      canContinue |= (jB[iA] != endB[iA]);
    }
 
    while (canContinue) {
      // Calculate the minimal indexB: selection of a candidate faceB
      Integer minIndexB = voidIndex;
      for (Integer iA = 0; iA < split_sizeA; ++iA) {
        Integer indexB = voidIndex;
        if (jB[iA] != endB[iA])
          indexB = jB[iA]->indexB();
        if (indexB != voidIndex) {
          if (minIndexB == voidIndex)
            minIndexB = indexB;
          else
            minIndexB = math::min(minIndexB, indexB);
        }
      }
 
      // Operator which convert Iter to Face (according to indexVoid special case)
      struct ConvertFace
      {
        ConvertFace(const Array<Face>& face_array, const Integer voidIndex)
        : m_face_array(face_array)
        , m_void_index(voidIndex)
        {}
        Face operator()(Integer index) const
        {
          if (index == m_void_index)
            return Face();
          else
            return m_face_array[index];
        }
 
       private:
 
        const Array<Face>& m_face_array;
        const Integer m_void_index;
      } convertFace(face_arrayB, voidIndex);
 
      // Use this minIndexB to get faceB
      Face faceB = convertFace(minIndexB);
 
      // Current accumulator for co-refinement values for faceA/faceB
      FaceFaceContact c(faceA, faceB);
 
      // Collection of information associated with the current faceB
      canContinue = false; // recalculate possible end
      for (Integer iA = 0; iA < split_sizeA; ++iA) {
        while (jB[iA] != endB[iA] and convertFace(jB[iA]->indexB()) == faceB) {
          const Handle& handle = *jB[iA];
          Real3 normal = convertGKVector(handle.normalA());
          Real3 center = convertGKVector(handle.centreA());
          c.centerA += center * math::normeR3(normal);
          if (face_reorientA[indexA + iA])
            c.normalA -= normal;
          else
            c.normalA += normal;
 
          if (not faceB.null()) {
            Real3 normal = convertGKVector(handle.normalB());
            Real3 center = convertGKVector(handle.centreB());
            c.centerB += center * math::normeR3(normal);
            if (face_reorientB[jB[iA]->indexB()])
              c.normalB -= normal;
            else
              c.normalB += normal;
          }
          ++jB[iA];
        }
        canContinue |= (jB[iA] != endB[iA]);
      }
      c.centerA /= math::normeR3(c.normalA);
      if (not faceB.null())
        c.centerB /= math::normeR3(c.normalB);
      contact.add(c);
      ++count;
    }
    indexA += split_sizeA;
  }
 
  { // case for indexA == voidIndex
    Iter jB = s.beginB(voidIndex);
    Iter endB = s.endB(voidIndex);
    while (jB != endB) {
      ARCANE_ASSERT((jB->indexB() != voidIndex), ("By construction indexB cannot be voidIndex when indexA is voidIndex"));
      Face faceB = face_arrayB[jB->indexB()];
 
      // Current accumulator for co-refinement values for faceA/faceB
      FaceFaceContact c(Face(), faceB);
 
      // Collection of information associated with the current faceB
      while (jB != endB and face_arrayB[jB->indexB()] == faceB) {
        const Handle& handle = *jB;
        Real3 normal = convertGKVector(handle.normalB());
        Real3 center = convertGKVector(handle.centreB());
        c.centerB += center * math::normeR3(normal);
        if (face_reorientB[jB->indexB()])
          c.normalB -= normal;
        else
          c.normalB += normal;
        ++jB;
      }
      c.centerB /= math::normeR3(c.normalB);
      contact.add(c);
      ++count;
    }
  }
 
  info() << "Co-refinement"
         << " : maxDeviation = " << s.maxDeviation()
         << " ; meanDeviation = " << s.meanDeviation()
         << " ; connections = " << count;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
GeometryKernelSurfaceImpl*
GeometryKernelSurfaceToolsService::
castSurface(ISurface* s) const
{
  GeometryKernelSurfaceImpl* gks = dynamic_cast<GeometryKernelSurfaceImpl*>(s);
  if (gks == NULL)
    fatal() << "Cannot use non own ISurface implementation";
  return gks;
}
 
/*---------------------------------------------------------------------------*/
 
void GeometryKernelSurfaceToolsService::
buildFaceGroupSurface(FaceGroup group, GeometryKernelSurfaceImpl* surface_impl) const
{
  ::buildFaceGroupSurface(group,
                          surface_impl->m_triangulation,
                          surface_impl->m_node_array,
                          surface_impl->m_face_array,
                          surface_impl->m_face_reorient,
                          surface_impl->m_mean_normal);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
ARCANE_REGISTER_SERVICE_GEOMETRYKERNELSURFACETOOLS(GeometryKernelSurfaceTools, GeometryKernelSurfaceToolsService);