]> SALOME platform Git repositories - tools/medcoupling.git/commitdiff
Salome HOME
Implementation of MEDCouplingPointSet.computeDiameterField (EDF10718) agy/ComputeDiameter
authorAnthony Geay <anthony.geay@edf.fr>
Wed, 29 Apr 2015 06:58:30 +0000 (08:58 +0200)
committerAnthony Geay <anthony.geay@edf.fr>
Wed, 29 Apr 2015 06:58:30 +0000 (08:58 +0200)
12 files changed:
src/INTERP_KERNEL/CMakeLists.txt
src/INTERP_KERNEL/CellModel.cxx
src/INTERP_KERNEL/CellModel.hxx
src/INTERP_KERNEL/DiameterCalculator.cxx [new file with mode: 0644]
src/INTERP_KERNEL/DiameterCalculator.hxx [new file with mode: 0644]
src/MEDCoupling/MEDCoupling1GTUMesh.cxx
src/MEDCoupling/MEDCoupling1GTUMesh.hxx
src/MEDCoupling/MEDCouplingPointSet.hxx
src/MEDCoupling/MEDCouplingUMesh.cxx
src/MEDCoupling/MEDCouplingUMesh.hxx
src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
src/MEDCoupling_Swig/MEDCouplingCommon.i

index 3523c9b85c680c6b8e4c3238121e171fdd7020d6..1fe07853bfbf2dc253a162d4732fe6e3765ae7b9 100644 (file)
@@ -26,6 +26,7 @@ SET(interpkernel_SOURCES
   TranslationRotationMatrix.cxx
   TetraAffineTransform.cxx
   CellModel.cxx
+  DiameterCalculator.cxx
   UnitTetraIntersectionBary.cxx
   InterpolationOptions.cxx
   BoxSplittingOptions.cxx
index 959162b0aed554ed5e8e212388780440d3a60684..7c01d990c40950c77fb4b5475098556de82f9964 100644 (file)
@@ -21,6 +21,7 @@
 #include "CellModel.hxx"
 
 #include "InterpKernelException.hxx"
+#include "DiameterCalculator.hxx"
 
 #include <algorithm>
 #include <sstream>
@@ -735,5 +736,67 @@ namespace INTERP_KERNEL
           }
       }
   }
-
+  
+  DiameterCalculator *CellModel::buildInstanceOfDiameterCalulator(int spaceDim) const
+  {
+    switch(_type)
+      {
+      case NORM_TRI3:
+        {
+          switch(spaceDim)
+            {
+            case 2:
+              return new DiameterCalulatorTRI3S2;
+            case 3:
+              return new DiameterCalulatorTRI3S3;
+            default:
+              throw Exception("CellModel::buildInstanceOfDiameterCalulator : For TRI3 only space dimension 2 and 3 implemented !");
+            }
+          break;
+        }
+      case NORM_QUAD4:
+        {
+          switch(spaceDim)
+            {
+            case 2:
+              return new DiameterCalulatorQUAD4S2;
+            case 3:
+              return new DiameterCalulatorQUAD4S3;
+            default:
+              throw Exception("CellModel::buildInstanceOfDiameterCalulator : For QUAD4 only space dimension 2 and 3 implemented !");
+            }
+          break;
+        }
+      case NORM_TETRA4:
+        {
+          if(spaceDim==3)
+            return new DiameterCalulatorTETRA4;
+          else
+            throw Exception("CellModel::buildInstanceOfDiameterCalulator : For TETRA4 space dimension 3 expected !");
+        }
+      case NORM_HEXA8:
+        {
+          if(spaceDim==3)
+            return new DiameterCalulatorHEXA8;
+          else
+            throw Exception("CellModel::buildInstanceOfDiameterCalulator : For HEXA8 space dimension 3 expected !");
+        }
+      case NORM_PENTA6:
+        {
+          if(spaceDim==3)
+            return new DiameterCalulatorPENTA6;
+          else
+            throw Exception("CellModel::buildInstanceOfDiameterCalulator : For PENTA6 space dimension 3 expected !");
+        }
+      case NORM_PYRA5:
+        {
+          if(spaceDim==3)
+            return new DiameterCalulatorPYRA5;
+          else
+            throw Exception("CellModel::buildInstanceOfDiameterCalulator : For PYRA5 space dimension 3 expected !");
+        }
+      default:
+        throw Exception("CellModel::buildInstanceOfDiameterCalulator : implemented only for TRI3, QUAD4, TETRA4, HEXA8, PENTA6, PYRA5 !");
+      }
+  }
 }
index 6b709012efe10f22e9c0861626f3754a6674c6ea..d1b5e1d956657b36677a3d51bc0304de40d9759c 100644 (file)
@@ -29,6 +29,8 @@
 
 namespace INTERP_KERNEL
 {
+  class DiameterCalculator;
+  
   /*!
    * This class descibes all static elements (different from polygons and polyhedron) 3D, 2D and 1D.
    */
@@ -74,6 +76,7 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT unsigned fillSonEdgesNodalConnectivity3D(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, NormalizedCellType& typeOfSon) const;
     INTERPKERNEL_EXPORT void changeOrientationOf2D(int *nodalConn, unsigned int sz) const;
     INTERPKERNEL_EXPORT void changeOrientationOf1D(int *nodalConn, unsigned int sz) const;
+    INTERPKERNEL_EXPORT DiameterCalculator *buildInstanceOfDiameterCalulator(int spaceDim) const;
   private:
     bool _dyn;
     bool _quadratic;
diff --git a/src/INTERP_KERNEL/DiameterCalculator.cxx b/src/INTERP_KERNEL/DiameterCalculator.cxx
new file mode 100644 (file)
index 0000000..a273959
--- /dev/null
@@ -0,0 +1,378 @@
+// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+//
+// 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.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// Author : Anthony Geay (EDF R&D)
+
+#include "DiameterCalculator.hxx"
+#include "InterpKernelException.hxx"
+#include "CellModel.hxx"
+
+#include <algorithm>
+#include <sstream>
+#include <cmath>
+
+using namespace INTERP_KERNEL;
+
+NormalizedCellType DiameterCalulatorTRI3S2::TYPE=NORM_TRI3;
+
+NormalizedCellType DiameterCalulatorTRI3S3::TYPE=NORM_TRI3;
+
+NormalizedCellType DiameterCalulatorQUAD4S2::TYPE=NORM_QUAD4;
+
+NormalizedCellType DiameterCalulatorQUAD4S3::TYPE=NORM_QUAD4;
+
+NormalizedCellType DiameterCalulatorTETRA4::TYPE=NORM_TETRA4;
+
+NormalizedCellType DiameterCalulatorHEXA8::TYPE=NORM_HEXA8;
+
+NormalizedCellType DiameterCalulatorPENTA6::TYPE=NORM_PENTA6;
+
+NormalizedCellType DiameterCalulatorPYRA5::TYPE=NORM_PYRA5;
+
+inline double SqNormV2(const double tab[2])
+{
+  return tab[0]*tab[0]+tab[1]*tab[1];
+}
+
+inline void DiffV2(const double a[2], const double b[2], double c[2])
+{
+  c[0]=a[0]-b[0];
+  c[1]=a[1]-b[1];
+}
+
+inline double SqNormV3(const double tab[3])
+{
+  return tab[0]*tab[0]+tab[1]*tab[1]+tab[2]*tab[2];
+}
+
+inline void DiffV3(const double a[3], const double b[3], double c[3])
+{
+  c[0]=a[0]-b[0];
+  c[1]=a[1]-b[1];
+  c[2]=a[2]-b[2];
+}
+
+template<class Evaluator>
+void ComputeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr)
+{
+  Evaluator evtor;
+  NormalizedCellType ct(Evaluator::TYPE);
+  int cti((int) ct);
+  for(const int *it=bgIds;it!=endIds;it++)
+    {
+      int offset(indPtr[*it]);
+      if(connPtr[offset]==cti)
+        resPtr[*it]=evtor.computeForOneCellInternal(connPtr+offset+1,connPtr+indPtr[(*it)+1],coordsPtr);
+      else
+        {
+          std::ostringstream oss; oss << "DiameterCalculator::computeForListOfCellIdsUMeshFrmt : invalid nodal connectivity format at cell # " << *it <<  " !";
+          throw Exception(oss.str().c_str());
+        }
+    }
+}
+
+template<class Evaluator>
+void ComputeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr)
+{
+  Evaluator evtor;
+  NormalizedCellType ct(Evaluator::TYPE);
+  int cti((int) ct);
+  for(int it=bgId;it<endId;it++)
+    {
+      int offset(indPtr[it]);
+      if(connPtr[offset]==cti)
+        resPtr[it]=evtor.computeForOneCellInternal(connPtr+offset+1,connPtr+indPtr[it+1],coordsPtr);
+      else
+        {
+          std::ostringstream oss; oss << "DiameterCalculator::computeForListOfCellIdsUMeshFrmt : invalid nodal connectivity format at cell # " << it <<  " !";
+          throw Exception(oss.str().c_str());
+        }
+    }
+}
+
+template<class Evaluator>
+void ComputeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr)
+{
+  Evaluator evtor;
+  NormalizedCellType ct(Evaluator::TYPE);
+  const CellModel& cm(CellModel::GetCellModel(ct));
+  unsigned nbNodes(cm.getNumberOfNodes());
+  const int *ptr(connPtr);
+  for(int i=0;i<nbOfCells;i++,ptr+=nbNodes,resPtr++)
+    *resPtr=evtor.computeForOneCellInternal(ptr,ptr+nbNodes,coordsPtr);
+}
+
+double DiameterCalulatorTRI3S2::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==3)
+    {
+      const double *a(coordsPtr+2*bg[0]),*b(coordsPtr+2*bg[1]),*c(coordsPtr+2*bg[2]);
+      double l0[2],l1[2],l2[2];
+      DiffV2(a,b,l0);
+      DiffV2(a,c,l1);
+      DiffV2(b,c,l2);
+      double res(std::max(SqNormV2(l0),SqNormV2(l1)));
+      res=std::max(res,SqNormV2(l2));
+      return std::sqrt(res);
+    }
+  else
+    throw Exception("DiameterCalulatorTRI3S2::computeForOneCellInternal : input connectivity must be of size 3 !");
+}
+
+void DiameterCalulatorTRI3S2::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S2>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTRI3S2::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S2>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTRI3S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorTRI3S2>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorTRI3S3::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==3)
+    {
+      const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]);
+      double l0[3],l1[3],l2[3];
+      DiffV3(a,b,l0);
+      DiffV3(a,c,l1);
+      DiffV3(b,c,l2);
+      double res(std::max(SqNormV3(l0),SqNormV3(l1)));
+      res=std::max(res,SqNormV3(l2));
+      return std::sqrt(res);
+    }
+  else
+    throw Exception("DiameterCalulatorTRI3S2::computeForOneCellInternal : input connectivity must be of size 3 !");
+}
+
+void DiameterCalulatorTRI3S3::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S3>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTRI3S3::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTRI3S3>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTRI3S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorTRI3S3>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorQUAD4S2::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==4)
+    {
+      const double *a(coordsPtr+2*bg[0]),*b(coordsPtr+2*bg[1]),*c(coordsPtr+2*bg[2]),*d(coordsPtr+2*bg[3]);
+      double l0[2],l1[2];
+      DiffV2(a,c,l0);
+      DiffV2(b,d,l1);
+      return std::sqrt(std::max(SqNormV2(l0),SqNormV2(l1)));
+    }
+  else
+    throw Exception("DiameterCalulatorQUAD4S2::computeForOneCellInternal : input connectivity must be of size 4 !");
+}
+
+void DiameterCalulatorQUAD4S2::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S2>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorQUAD4S2::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S2>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorQUAD4S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorQUAD4S2>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorQUAD4S3::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==4)
+    {
+      const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]),*d(coordsPtr+3*bg[3]);
+      double l0[3],l1[3];
+      DiffV3(a,c,l0);
+      DiffV3(b,d,l1);
+      return std::sqrt(std::max(SqNormV3(l0),SqNormV3(l1)));
+    }
+  else
+    throw Exception("DiameterCalulatorQUAD4S3::computeForOneCellInternal : input connectivity must be of size 4 !");
+}
+
+void DiameterCalulatorQUAD4S3::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S3>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorQUAD4S3::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorQUAD4S3>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorQUAD4S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorQUAD4S3>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorTETRA4::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==4)
+    {
+      const double *a(coordsPtr+3*bg[0]),*b(coordsPtr+3*bg[1]),*c(coordsPtr+3*bg[2]),*d(coordsPtr+3*bg[3]);
+      double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
+      DiffV3(a,b,l0);
+      DiffV3(a,c,l1);
+      DiffV3(b,c,l2);
+      DiffV3(a,d,l3);
+      DiffV3(b,d,l4);
+      DiffV3(c,d,l5);
+      double tmp[6];
+      tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
+      return std::sqrt(*std::max_element(tmp,tmp+6));
+    }
+  else
+    throw Exception("DiameterCalulatorTETRA4::computeForOneCellInternal : input connectivity must be of size 4 !");
+}
+
+void DiameterCalulatorTETRA4::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorTETRA4>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTETRA4::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorTETRA4>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorTETRA4::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorTETRA4>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorHEXA8::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==8)
+    {
+      const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]),*p5(coordsPtr+3*bg[5]),*p6(coordsPtr+3*bg[6]),*p7(coordsPtr+3*bg[7]);
+      double l0[3],l1[3],l2[3],l3[3];
+      DiffV3(p0,p6,l0);
+      DiffV3(p1,p7,l1);
+      DiffV3(p2,p4,l2);
+      DiffV3(p3,p5,l3);
+      double tmp[4];
+      tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3);
+      return std::sqrt(*std::max_element(tmp,tmp+4));
+    }
+  else
+    throw Exception("DiameterCalulatorHEXA8::computeForOneCellInternal : input connectivity must be of size 8 !");
+}
+
+void DiameterCalulatorHEXA8::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorHEXA8>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorHEXA8::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorHEXA8>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorHEXA8::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorHEXA8>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorPENTA6::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==6)
+    {
+      const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]),*p5(coordsPtr+3*bg[5]);
+      double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
+      DiffV3(p0,p4,l0);
+      DiffV3(p1,p3,l1);
+      DiffV3(p1,p5,l2);
+      DiffV3(p2,p4,l3);
+      DiffV3(p0,p5,l4);
+      DiffV3(p2,p3,l5);
+      double tmp[6];
+      tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
+      return std::sqrt(*std::max_element(tmp,tmp+6));
+    }
+  else
+    throw Exception("DiameterCalulatorPENTA6::computeForOneCellInternal : input connectivity must be of size 6 !");
+}
+
+void DiameterCalulatorPENTA6::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorPENTA6>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorPENTA6::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorPENTA6>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorPENTA6::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorPENTA6>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
+
+double DiameterCalulatorPYRA5::computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const
+{
+  if(std::distance(bg,endd)==5)
+    {
+      const double *p0(coordsPtr+3*bg[0]),*p1(coordsPtr+3*bg[1]),*p2(coordsPtr+3*bg[2]),*p3(coordsPtr+3*bg[3]),*p4(coordsPtr+3*bg[4]);
+      double l0[3],l1[3],l2[3],l3[3],l4[3],l5[3];
+      DiffV3(p0,p2,l0);
+      DiffV3(p1,p3,l1);
+      DiffV3(p0,p4,l2);
+      DiffV3(p1,p4,l3);
+      DiffV3(p2,p4,l4);
+      DiffV3(p3,p4,l5);
+      double tmp[6];
+      tmp[0]=SqNormV3(l0); tmp[1]=SqNormV3(l1); tmp[2]=SqNormV3(l2); tmp[3]=SqNormV3(l3); tmp[4]=SqNormV3(l4); tmp[5]=SqNormV3(l5);
+      return std::sqrt(*std::max_element(tmp,tmp+6));
+    }
+  else
+    throw Exception("DiameterCalulatorPYRA5::computeForOneCellInternal : input connectivity must be of size 5 !");
+}
+
+void DiameterCalulatorPYRA5::computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForListOfCellIdsUMeshFrmt<DiameterCalulatorPYRA5>(bgIds,endIds,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorPYRA5::computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeForRangeOfCellIdsUMeshFrmt<DiameterCalulatorPYRA5>(bgId,endId,indPtr,connPtr,coordsPtr,resPtr);
+}
+
+void DiameterCalulatorPYRA5::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const
+{
+  ComputeFor1SGTUMeshFrmt<DiameterCalulatorPYRA5>(nbOfCells,connPtr,coordsPtr,resPtr);
+}
diff --git a/src/INTERP_KERNEL/DiameterCalculator.hxx b/src/INTERP_KERNEL/DiameterCalculator.hxx
new file mode 100644 (file)
index 0000000..865e866
--- /dev/null
@@ -0,0 +1,146 @@
+// Copyright (C) 2007-2015  CEA/DEN, EDF R&D
+//
+// 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.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// Author : Anthony Geay (EDF R&D)
+
+#ifndef __DIAMETERCALCULATOR_HXX__
+#define __DIAMETERCALCULATOR_HXX__
+
+#include "INTERPKERNELDefines.hxx"
+
+#include "NormalizedGeometricTypes"
+
+namespace INTERP_KERNEL
+{
+  class DiameterCalculator
+  {
+  public:
+    INTERPKERNEL_EXPORT virtual ~DiameterCalculator() { }
+    INTERPKERNEL_EXPORT virtual NormalizedCellType getType() const = 0;
+    INTERPKERNEL_EXPORT virtual double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const = 0;
+    INTERPKERNEL_EXPORT virtual void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const = 0;
+    INTERPKERNEL_EXPORT virtual void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const = 0;
+    INTERPKERNEL_EXPORT virtual void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const = 0;
+  };
+
+  class DiameterCalulatorTRI3S2 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorTRI3S3 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorQUAD4S2 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorQUAD4S3 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorTETRA4 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorHEXA8 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorPENTA6 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+
+  class DiameterCalulatorPYRA5 : public DiameterCalculator
+  {
+  public:
+    NormalizedCellType getType() const { return TYPE; }
+    double computeForOneCell(const int *bg, const int *endd, const double *coordsPtr) const { return computeForOneCellInternal(bg,endd,coordsPtr); }
+    double computeForOneCellInternal(const int *bg, const int *endd, const double *coordsPtr) const;
+    void computeForListOfCellIdsUMeshFrmt(const int *bgIds, const int *endIds, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeForRangeOfCellIdsUMeshFrmt(int bgId, int endId, const int *indPtr, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+    void computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const;
+  public:
+    static NormalizedCellType TYPE;
+  };
+}
+
+#endif
index 343076dfd29a143f7a870ec64d7ac8fb8f1ff001..eb0db7da861b956ef7161534d46e55257b16ce58 100644 (file)
@@ -24,6 +24,8 @@
 #include "MEDCouplingCMesh.hxx"
 
 #include "SplitterTetra.hxx"
+#include "DiameterCalculator.hxx"
+#include "InterpKernelAutoPtr.hxx"
 
 using namespace ParaMEDMEM;
 
@@ -2113,6 +2115,26 @@ DataArrayDouble *MEDCoupling1SGTUMesh::getBoundingBoxForBBTree(double arcDetEps)
   return ret.retn();
 }
 
+/*!
+ * Returns the cell field giving for each cell in \a this its diameter. Diameter means the max length of all possible SEG2 in the cell.
+ *
+ * \return a new instance of field containing the result. The returned instance has to be deallocated by the caller.
+ */
+MEDCouplingFieldDouble *MEDCoupling1SGTUMesh::computeDiameterField() const
+{
+  checkFullyDefined();
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
+  int nbCells(getNumberOfCells());
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::New());
+  arr->alloc(nbCells,1);
+  INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::DiameterCalculator> dc(_cm->buildInstanceOfDiameterCalulator(getSpaceDimension()));
+  dc->computeFor1SGTUMeshFrmt(nbCells,_conn->begin(),getCoords()->begin(),arr->getPointer());
+  ret->setMesh(this);
+  ret->setArray(arr);
+  ret->setName("Diameter");
+  return ret.retn();
+}
+
 //== 
 
 MEDCoupling1DGTUMesh *MEDCoupling1DGTUMesh::New()
@@ -3594,6 +3616,16 @@ DataArrayDouble *MEDCoupling1DGTUMesh::getBoundingBoxForBBTree(double arcDetEps)
   return ret.retn();
 }
 
+/*!
+ * Returns the cell field giving for each cell in \a this its diameter. Diameter means the max length of all possible SEG2 in the cell.
+ *
+ * \return a new instance of field containing the result. The returned instance has to be deallocated by the caller.
+ */
+MEDCouplingFieldDouble *MEDCoupling1DGTUMesh::computeDiameterField() const
+{
+  throw INTERP_KERNEL::Exception("MEDCoupling1DGTUMesh::computeDiameterField : not implemented yet for dynamic types !");
+}
+
 std::vector<int> MEDCoupling1DGTUMesh::BuildAPolygonFromParts(const std::vector< std::vector<int> >& parts)
 {
   std::vector<int> ret;
index 3ccb85cb51ade26f3cec4c832ee99d8040bf2edc..9928bc48f23d3681d09154c5b1e63b71fe43aa25 100644 (file)
@@ -136,6 +136,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const;
     MEDCOUPLING_EXPORT int getNumberOfNodesInCell(int cellId) const;
     MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *computeDiameterField() const;
     // overload of MEDCoupling1GTUMesh
     MEDCOUPLING_EXPORT void checkCoherencyOfConnectivity() const;
     MEDCOUPLING_EXPORT void allocateCells(int nbOfCells=0);
@@ -229,6 +230,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, DataArrayInt *&cellIdsKeptArr) const;
     MEDCOUPLING_EXPORT int getNumberOfNodesInCell(int cellId) const;
     MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *computeDiameterField() const;
     // overload of MEDCoupling1GTUMesh
     MEDCOUPLING_EXPORT void checkCoherencyOfConnectivity() const;
     MEDCOUPLING_EXPORT void allocateCells(int nbOfCells=0);
index 08281ec962bf71b9b9407246b130f8aa545f4aff..21c830d2f2abde8a1c6f75adf38f1fe42316c64e 100644 (file)
@@ -139,6 +139,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT virtual DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *getCellsInBoundingBox(const double *bbox, double eps) const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps) = 0;
+    MEDCOUPLING_EXPORT virtual MEDCouplingFieldDouble *computeDiameterField() const = 0;
     MEDCOUPLING_EXPORT virtual DataArrayInt *zipCoordsTraducer();
     MEDCOUPLING_EXPORT virtual DataArrayInt *zipConnectivityTraducer(int compType, int startCellId=0);
     MEDCOUPLING_EXPORT virtual bool areAllNodesFetched() const;
index f213a9d42d53de8b084d6e3fb19c1913bfbc112e..b8421f25122708abeb94274c1d0dd04044e367c5 100644 (file)
@@ -29,6 +29,7 @@
 #include "BBTree.txx"
 #include "BBTreeDst.txx"
 #include "SplitterTetra.hxx"
+#include "DiameterCalculator.hxx"
 #include "DirectedBoundingBox.hxx"
 #include "InterpKernelMatrixTools.hxx"
 #include "InterpKernelMeshQuality.hxx"
@@ -3145,16 +3146,7 @@ MEDCouplingUMesh::~MEDCouplingUMesh()
  */
 void MEDCouplingUMesh::computeTypes()
 {
-  if(_nodal_connec && _nodal_connec_index)
-    {
-      _types.clear();
-      const int *conn=_nodal_connec->getConstPointer();
-      const int *connIndex=_nodal_connec_index->getConstPointer();
-      int nbOfElem=_nodal_connec_index->getNbOfElems()-1;
-      if (nbOfElem > 0)
-        for(const int *pt=connIndex;pt !=connIndex+nbOfElem;pt++)
-          _types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
-    }
+  ComputeAllTypesInternal(_types,_nodal_connec,_nodal_connec_index);
 }
 
 /*!
@@ -6515,6 +6507,34 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const
   return ret.retn();
 }
 
+/*!
+ * Returns the cell field giving for each cell in \a this its diameter. Diameter means the max length of all possible SEG2 in the cell.
+ *
+ * \return a new instance of field containing the result. The returned instance has to be deallocated by the caller.
+ *
+ * \sa getSkewField, getWarpField, getAspectRatioField, getEdgeRatioField
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::computeDiameterField() const
+{
+  checkCoherency();
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
+  ret->setMesh(this);
+  std::set<INTERP_KERNEL::NormalizedCellType> types;
+  ComputeAllTypesInternal(types,_nodal_connec,_nodal_connec_index);
+  int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells());
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::New());
+  arr->alloc(nbCells,1);
+  for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++)
+    {
+      INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::DiameterCalculator> dc(INTERP_KERNEL::CellModel::GetCellModel(*it).buildInstanceOfDiameterCalulator(spaceDim));
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cellIds(giveCellsWithType(*it));
+      dc->computeForListOfCellIdsUMeshFrmt(cellIds->begin(),cellIds->end(),_nodal_connec_index->begin(),_nodal_connec->begin(),getCoords()->begin(),arr->getPointer());
+    }
+  ret->setArray(arr);
+  ret->setName("Diameter");
+  return ret.retn();
+}
+
 /*!
  * This method aggregate the bbox of each cell and put it into bbox parameter.
  * 
@@ -11535,6 +11555,19 @@ int MEDCouplingUMesh::split2DCellsQuadratic(const DataArrayInt *desc, const Data
   return addCoo->getNumberOfTuples();
 }
 
+void MEDCouplingUMesh::ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex)
+{
+  if(nodalConnec && nodalConnecIndex)
+    {
+      types.clear();
+      const int *conn(nodalConnec->getConstPointer()),*connIndex(nodalConnecIndex->getConstPointer());
+      int nbOfElem(nodalConnecIndex->getNbOfElems()-1);
+      if(nbOfElem>0)
+        for(const int *pt=connIndex;pt!=connIndex+nbOfElem;pt++)
+          types.insert((INTERP_KERNEL::NormalizedCellType)conn[*pt]);
+    }
+}
+
 MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)),
     _own_cell(true),_cell_id(-1),_nb_cell(0)
 {
index a3972a295942aa6634fde4f2d2e2d53d0ee9d989..fff81f7c505aef62d06b4ac31b6a76d77b1817f7 100644 (file)
@@ -194,6 +194,7 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getAspectRatioField() const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getWarpField() const;
     MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getSkewField() const;
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *computeDiameterField() const;
     //utilities for MED File RW
     MEDCOUPLING_EXPORT std::vector<int> getDistributionOfTypes() const;
     MEDCOUPLING_EXPORT DataArrayInt *checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const;
@@ -333,6 +334,7 @@ namespace ParaMEDMEM
     void split2DCellsLinear(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI);
     int split2DCellsQuadratic(const DataArrayInt *desc, const DataArrayInt *descI, const DataArrayInt *subNodesInSeg, const DataArrayInt *subNodesInSegI, const DataArrayInt *mid, const DataArrayInt *midI);
     static bool Colinearize2DCell(const double *coords, const int *connBg, const int *connEnd, int offset, DataArrayInt *newConnOfCell, DataArrayDouble *appendedCoords);
+    static void ComputeAllTypesInternal(std::set<INTERP_KERNEL::NormalizedCellType>& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex);
   public:
     MEDCOUPLING_EXPORT static DataArrayInt *ComputeRangesFromTypeDistribution(const std::vector<int>& code);
     MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=24;
index c6bc0d7d085687dc7978453e9ae27eca387c7171..bcfa2c8635baedf4645ee872124f29c165f36a85 100644 (file)
@@ -16429,6 +16429,123 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         m2.zipCoords()
         self.assertTrue(m2.areAllNodesFetched())
         pass
+
+    def testMEDCouplingPointSetComputeDiameterField1(self):
+        arrX=DataArrayDouble([0.,1.1,1.7,2.1])
+        arrY=DataArrayDouble([0.,0.7,0.8,1.9])
+        arrZ=DataArrayDouble([0.,1.3,2.1,2.4])
+        m=MEDCouplingCMesh() ; m.setCoords(arrX,arrY,arrZ) ; m=m.buildUnstructured()
+        f=m.computeDiameterField()
+        f.checkCoherency()
+        exp=DataArrayDouble([1.8411952639521971,1.5937377450509227,1.5297058540778357,1.705872210923198,1.4352700094407325,1.3638181696985856,2.0273134932713295,1.8055470085267789,1.7492855684535902,1.5297058540778357,1.2206555615733703,1.1357816691600546,1.3638181696985856,1.004987562112089,0.9,1.7492855684535902,1.4866068747318506,1.4177446878757824,1.3379088160259651,0.9695359714832656,0.8602325267042626,1.1445523142259597,0.6782329983125266,0.5099019513592785,1.5842979517754858,1.2884098726725124,1.208304597359457])
+        self.assertTrue(exp.isEqual(f.getArray(),1e-12))
+        m1=m[::2]
+        m2=m[1::2]
+        m2.simplexize(PLANAR_FACE_5)
+        m3=MEDCouplingUMesh.MergeUMeshesOnSameCoords(m1,m2)
+        f=m3.computeDiameterField()
+        f.checkCoherency()
+        exp2=DataArrayDouble([1.8411952639521971,1.5297058540778357,1.4352700094407325,2.0273134932713295,1.7492855684535902,1.2206555615733703,1.3638181696985856,0.9,1.4866068747318506,1.3379088160259651,0.8602325267042626,0.6782329983125266,1.5842979517754858,1.208304597359457,1.47648230602334,1.47648230602334,1.47648230602334,1.47648230602334,1.47648230602334,1.7029386365926402,1.7029386365926402,1.7029386365926402,1.7029386365926402,1.7029386365926402,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.70293863659264,1.70293863659264,1.70293863659264,1.70293863659264,1.70293863659264,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.3601470508735445,1.063014581273465,1.063014581273465,1.063014581273465,1.063014581273465,1.063014581273465,1.0,1.0,1.0,1.0,1.0,1.5556349186104046,1.5556349186104046,1.5556349186104046,1.5556349186104046,1.5556349186104046,1.3601470508735443,1.3601470508735443,1.3601470508735443,1.3601470508735443,1.3601470508735443,0.9219544457292886,0.9219544457292886,0.9219544457292886,0.9219544457292886,0.9219544457292886,1.140175425099138,1.140175425099138,1.140175425099138,1.140175425099138,1.140175425099138,0.5,0.5,0.5,0.5,0.5,1.2529964086141667,1.2529964086141667,1.2529964086141667,1.2529964086141667,1.2529964086141667])
+        self.assertTrue(exp2.isEqual(f.getArray(),1e-12))
+        # TRI3 - spacedim = 2
+        coo=DataArrayDouble([(1,1),(5,1.9),(2.1,3)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_TRI3) ; m.setCoords(coo)
+        for c in [[0,1,2],[0,2,1],[2,1,0]]:
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],4.1,12)
+        # TRI3 - spacedim = 3
+        coo=DataArrayDouble([(1.3198537928820775,1.0991902391274959,-0.028645697595823361),(5.2486835106806335,2.2234012799688281,0.30368935050077939),(2.2973688139447361,3.1572023778066649,0.10937756365410012)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_TRI3) ; m.setCoords(coo)
+        for c in [[0,1,2],[0,2,1],[2,1,0]]:
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],4.1,12)
+        # QUAD4 - spacedim = 2
+        coo=DataArrayDouble([(0,2),(2,0),(6,4),(4,9)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_QUAD4) ; m.setCoords(coo)
+        exp3=sqrt(85.)
+        for delta in xrange(4):
+            c=[(elt+delta)%4 for elt in xrange(4)]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp3,12)
+            c.reverse()
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp3,12)
+        # QUAD4 - spacedim = 3
+        coo=DataArrayDouble([(0.26570992384234871,2.0405889913271817,-0.079134238105786903),(2.3739976619218064,0.15779148692781009,0.021842842914139737),(6.1207841448393197,4.3755532938679655,0.43666375769970678),(3.8363255342943359,9.2521096041694229,0.41551170895942313)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_QUAD4) ; m.setCoords(coo)
+        for delta in xrange(4):
+            c=[(elt+delta)%4 for elt in xrange(4)]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp3,12)
+            c.reverse()
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp3,12)
+        # PENTA6
+        # noise of coo=DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(0,0,2),(1,0,2),(0,1,2)]) + rotation([0.7,-1.2,0.6],[-4,-1,10],0.3)
+        coo=DataArrayDouble([(-0.28594726851554486,-0.23715005500928255,-0.10268080010083136),(0.6167364988633947,-0.008923258436324799,-0.08574087516687756),(-0.6132873463333834,0.6943403970881654,-0.2806118260037991),(-0.40705974936532896,-0.05868487929989308,1.7724055544436323),(0.5505955507861958,0.19145393798144705,1.8788156352163994),(-0.6092686217773406,0.812502961290914,1.685712743757831)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_PENTA6) ; m.setCoords(coo)
+        exp4=2.5041256256889888
+        self.assertAlmostEqual(exp4,coo.buildEuclidianDistanceDenseMatrix().getMaxValue()[0],12)# <- the definition of diameter
+        for delta in xrange(3):
+            c=[(elt+delta)%3 for elt in xrange(3)]
+            c+=[elt+3 for elt in c]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp4,12)
+            c.reverse()
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp4,12)
+        # HEXA8
+        # noise of coo=DataArrayDouble([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,2),(1,0,2),(1,1,2),(0,1,2)]) + rotation([0.7,-1.2,0.6],[-4,-1,10],0.3)
+        coo=DataArrayDouble([(-0.21266406388867243,-0.3049569460042527,-0.11012394815006032),(0.7641037943272584,-0.06990814759929553,-0.0909613877456491),(0.47406560768559974,0.8681310650341907,-0.2577311403703061),(-0.5136830410871793,0.644390554940524,-0.21319015989794698),(-0.4080167737381202,-0.12853761670628505,1.7869166291979348),(0.5650318811550441,0.20476257733110748,1.8140158890821603),(0.3230844436386215,1.1660778242678538,1.7175073141333406),(-0.6656588358432984,0.918357550969698,1.7566470691880265)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_HEXA8) ; m.setCoords(coo)
+        exp5=2.5366409441884215
+        self.assertAlmostEqual(exp5,coo.buildEuclidianDistanceDenseMatrix().getMaxValue()[0],12)# <- the definition of diameter
+        for delta in xrange(4):
+            c=[(elt+delta)%4 for elt in xrange(4)]
+            c+=[elt+4 for elt in c]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp5,12)
+            c.reverse()
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp5,12)
+        # PYRA5 (1) 5th node is further 
+        # noise of coo=DataArrayDouble([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0.5,0.5,2)]) + rotation([0.7,-1.2,0.6],[-4,-1,10],0.3)
+        coo=DataArrayDouble([(-0.31638393672228626,-0.3157865246451914,-0.12555467233075002),(0.7281379795666488,0.03836511217237115,-0.08431662762197323),(0.4757967840735147,0.8798897996143908,-0.2680890320119049),(-0.5386339871809047,0.5933159894201252,-0.2975311238319419),(0.012042592988768974,0.534282135495012,1.7859521682027926)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_PYRA5) ; m.setCoords(coo)
+        exp6=2.1558368027391386
+        self.assertAlmostEqual(exp6,coo.buildEuclidianDistanceDenseMatrix().getMaxValue()[0],12)# <- the definition of diameter
+        for delta in xrange(4):
+            c=[(elt+delta)%4 for elt in xrange(4)]
+            c+=[4]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp6,12)
+            pass
+        # PYRA5 (2) 5th node is closer
+        # noise of coo=DataArrayDouble([(0,0,0),(1,0,0),(1,1,0),(0,1,0),(0.5,0.5,0.1)]) + rotation([0.7,-1.2,0.6],[-4,-1,10],0.3)
+        coo=DataArrayDouble([(-0.31638393672228626,-0.3157865246451914,-0.12555467233075002),(0.7281379795666488,0.03836511217237115,-0.08431662762197323),(0.4757967840735147,0.8798897996143908,-0.2680890320119049),(-0.5386339871809047,0.5933159894201252,-0.2975311238319419),(0.092964408350795,0.33389670321297005,-0.10171764888060142)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_PYRA5) ; m.setCoords(coo)
+        exp7=1.4413563787228953
+        self.assertAlmostEqual(exp7,coo.buildEuclidianDistanceDenseMatrix().getMaxValue()[0],12)# <- the definition of diameter
+        for delta in xrange(4):
+            c=[(elt+delta)%4 for elt in xrange(4)]
+            c+=[4]
+            m.setNodalConnectivity(DataArrayInt(c))
+            self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp7,12)
+            pass
+        # TETRA4
+        # noise of coo=DataArrayDouble([(0,0,0),(1,0,0),(0,1,0),(1,1,1)]) + rotation([0.7,-1.2,0.6],[-4,-1,10],0.3)
+        coo=DataArrayDouble([(-0.2256894071281369,-0.27631691290428106,-0.20266086543995965),(0.655458695100186,-0.08173323565551605,-0.19254662462061933),(-0.49893490718947264,0.5848097154568599,-0.3039928255382145),(0.2988102920828487,1.0582266398878504,0.7347375047372364)])
+        m=MEDCoupling1SGTUMesh("mesh",NORM_TETRA4) ; m.setCoords(coo)
+        exp8=1.7131322579364157
+        self.assertAlmostEqual(exp8,coo.buildEuclidianDistanceDenseMatrix().getMaxValue()[0],12)# <- the definition of diameter
+        for c in [[0,1,2,3],[0,1,3,2],[0,3,2,1],[0,3,1,2]]:
+            for i in xrange(4):
+                m.setNodalConnectivity(DataArrayInt([(elt+i)%4 for elt in c]))
+                self.assertAlmostEqual(m.computeDiameterField().getArray()[0],exp8,12)
+                pass
+            pass
+        pass
+
     pass
 
 if __name__ == '__main__':
index 6926805aebff562030f0e6c9adf97c161ae4e017..22cd189214186fed689ec7b61c0878d2204cebb8 100644 (file)
@@ -260,6 +260,7 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingPointSet::getBoundingBoxForBBTree;
 %newobject ParaMEDMEM::MEDCouplingPointSet::computeFetchedNodeIds;
 %newobject ParaMEDMEM::MEDCouplingPointSet::ComputeNbOfInteractionsWithSrcCells;
+%newobject ParaMEDMEM::MEDCouplingPointSet::computeDiameterField;
 %newobject ParaMEDMEM::MEDCouplingPointSet::__getitem__;
 %newobject ParaMEDMEM::MEDCouplingUMesh::New;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getNodalConnectivity;
@@ -1174,6 +1175,7 @@ namespace ParaMEDMEM
       virtual DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
       virtual void renumberNodesWithOffsetInConn(int offset) throw(INTERP_KERNEL::Exception);
       virtual bool areAllNodesFetched() const throw(INTERP_KERNEL::Exception);
+      virtual MEDCouplingFieldDouble *computeDiameterField() const throw(INTERP_KERNEL::Exception);
       %extend 
          {
            std::string __str__() const throw(INTERP_KERNEL::Exception)