GeometryKernelSurfaceToolsService.cc

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2022 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
#include "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/IParallelMng.h>
#include <arcane/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 "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;
 
  // bidouille pour sauter les pb de pointeurs boost non voulus
  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 sur indexA est traité à part
  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 des positions courantes des sous-faces B associées à cette 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) {
      // Calcul de l'indexB minimal : selection d'une faceB candidate
      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);
 
      // Utilise de cet minIndexB donne le faceB
      Face faceB = convertFace(minIndexB);
 
      // Current accumulator for co-refinement values for faceA/faceB
      FaceFaceContact c(faceA,faceB);
 
      // Collecte des informations associées à la faceB courante
      canContinue = false; // recalcul de la fin possible
      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);
 
      // Collecte des informations associées à la faceB courante
      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);