From: ageay <ageay>
Date: Mon, 18 Oct 2010 07:43:53 +0000 (+0000)
Subject: *** empty log message ***
X-Git-Tag: V5_1_5~11
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=98db8972a3182e1ef4d849e9745d54d5682696bd;p=tools%2Fmedcoupling.git

*** empty log message ***
---

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index ce0461314..4da00a715 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -25,6 +25,7 @@
 #include "PointLocatorAlgos.txx"
 #include "BBTree.txx"
 #include "DirectedBoundingBox.hxx"
+#include "InterpKernelMeshQuality.hxx"
 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
 
 #include <sstream>
@@ -2286,6 +2287,222 @@ void MEDCouplingUMesh::getFastAveragePlaneOfThis(double *vec, double *pos) const
   std::copy(coordsPtr+3*conn[1],coordsPtr+3*conn[1]+3,pos);
 }
 
+/*!
+ * The returned newly created field has to be managed by the caller.
+ * This method returns a field on cell with no time lying on 'this'. The meshdimension and spacedimension of this are expected to be both in [2,3]. If not an exception will be thrown.
+ * This method for the moment only deals with NORM_TRI3, NORM_QUAD4 and NORM_TETRA4 geometric types.
+ * If a cell has an another type an exception will be thrown.
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::getEdgeRatioField() const throw(INTERP_KERNEL::Exception)
+{
+  checkCoherency();
+  int spaceDim=getSpaceDimension();
+  int meshDim=getMeshDimension();
+  if(spaceDim!=2 && spaceDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : SpaceDimension must be equal to 2 or 3 !");
+  if(meshDim!=2 && meshDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : MeshDimension must be equal to 2 or 3 !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  ret->setMesh(this);
+  int nbOfCells=getNumberOfCells();
+  DataArrayDouble *arr=DataArrayDouble::New();
+  arr->alloc(nbOfCells,1);
+  double *pt=arr->getPointer();
+  ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
+  arr->decrRef();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  const double *coo=_coords->getConstPointer();
+  double tmp[12];
+  for(int i=0;i<nbOfCells;i++,pt++)
+    {
+      INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
+      switch(t)
+        {
+          case INTERP_KERNEL::NORM_TRI3:
+            {
+              fillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::triEdgeRatio(tmp);
+              break;
+            }
+          case INTERP_KERNEL::NORM_QUAD4:
+            {
+              fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::quadEdgeRatio(tmp);
+              break;
+            }
+          case INTERP_KERNEL::NORM_TETRA4:
+            {
+              fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::tetraEdgeRatio(tmp);
+              break;
+            }
+        default:
+          throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getEdgeRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
+        }
+      conn+=connI[i+1]-connI[i];
+    }
+  ret->setName("EdgeRatio");
+  ret->incrRef();
+  return ret;
+}
+
+/*!
+ * The returned newly created field has to be managed by the caller.
+ * This method returns a field on cell with no time lying on 'this'. The meshdimension and spacedimension of this are expected to be both in [2,3]. If not an exception will be thrown.
+ * This method for the moment only deals with NORM_TRI3, NORM_QUAD4 and NORM_TETRA4 geometric types.
+ * If a cell has an another type an exception will be thrown.
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::getAspectRatioField() const throw(INTERP_KERNEL::Exception)
+{
+  checkCoherency();
+  int spaceDim=getSpaceDimension();
+  int meshDim=getMeshDimension();
+  if(spaceDim!=2 && spaceDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : SpaceDimension must be equal to 2 or 3 !");
+  if(meshDim!=2 && meshDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : MeshDimension must be equal to 2 or 3 !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  ret->setMesh(this);
+  int nbOfCells=getNumberOfCells();
+  DataArrayDouble *arr=DataArrayDouble::New();
+  arr->alloc(nbOfCells,1);
+  double *pt=arr->getPointer();
+  ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
+  arr->decrRef();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  const double *coo=_coords->getConstPointer();
+  double tmp[12];
+  for(int i=0;i<nbOfCells;i++,pt++)
+    {
+      INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
+      switch(t)
+        {
+          case INTERP_KERNEL::NORM_TRI3:
+            {
+              fillInCompact3DMode(spaceDim,3,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::triAspectRatio(tmp);
+              break;
+            }
+          case INTERP_KERNEL::NORM_QUAD4:
+            {
+              fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::quadAspectRatio(tmp);
+              break;
+            }
+          case INTERP_KERNEL::NORM_TETRA4:
+            {
+              fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::tetraAspectRatio(tmp);
+              break;
+            }
+        default:
+          throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getAspectRatioField : A cell with not manged type (NORM_TRI3, NORM_QUAD4 and NORM_TETRA4) has been detected !");
+        }
+      conn+=connI[i+1]-connI[i];
+    }
+  ret->setName("AspectRatio");
+  ret->incrRef();
+  return ret;
+}
+
+/*!
+ * The returned newly created field has to be managed by the caller.
+ * This method returns a field on cell with no time lying on 'this'. The meshdimension must be equal to 2 and the spacedimension must be equal to 3. If not an exception will be thrown.
+ * This method for the moment only deals with NORM_QUAD4 geometric type.
+ * If a cell has an another type an exception will be thrown.
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::getWarpField() const throw(INTERP_KERNEL::Exception)
+{
+  checkCoherency();
+  int spaceDim=getSpaceDimension();
+  int meshDim=getMeshDimension();
+  if(spaceDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : SpaceDimension must be equal to 3 !");
+  if(meshDim!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : MeshDimension must be equal to 2 !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  ret->setMesh(this);
+  int nbOfCells=getNumberOfCells();
+  DataArrayDouble *arr=DataArrayDouble::New();
+  arr->alloc(nbOfCells,1);
+  double *pt=arr->getPointer();
+  ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
+  arr->decrRef();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  const double *coo=_coords->getConstPointer();
+  double tmp[12];
+  for(int i=0;i<nbOfCells;i++,pt++)
+    {
+      INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
+      switch(t)
+        {
+          case INTERP_KERNEL::NORM_QUAD4:
+            {
+              fillInCompact3DMode(3,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::quadWarp(tmp);
+              break;
+            }
+        default:
+          throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getWarpField : A cell with not manged type (NORM_QUAD4) has been detected !");
+        }
+      conn+=connI[i+1]-connI[i];
+    }
+  ret->setName("Warp");
+  ret->incrRef();
+  return ret;
+}
+
+/*!
+ * The returned newly created field has to be managed by the caller.
+ * This method returns a field on cell with no time lying on 'this'. The meshdimension must be equal to 2 and the spacedimension must be equal to 3. If not an exception will be thrown.
+ * This method for the moment only deals with NORM_QUAD4 geometric type.
+ * If a cell has an another type an exception will be thrown.
+ */
+MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const throw(INTERP_KERNEL::Exception)
+{
+  checkCoherency();
+  int spaceDim=getSpaceDimension();
+  int meshDim=getMeshDimension();
+  if(spaceDim!=3)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : SpaceDimension must be equal to 3 !");
+  if(meshDim!=2)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : MeshDimension must be equal to 2 !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  ret->setMesh(this);
+  int nbOfCells=getNumberOfCells();
+  DataArrayDouble *arr=DataArrayDouble::New();
+  arr->alloc(nbOfCells,1);
+  double *pt=arr->getPointer();
+  ret->setArray(arr);//In case of throw to avoid mem leaks arr will be used after decrRef.
+  arr->decrRef();
+  const int *conn=_nodal_connec->getConstPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  const double *coo=_coords->getConstPointer();
+  double tmp[12];
+  for(int i=0;i<nbOfCells;i++,pt++)
+    {
+      INTERP_KERNEL::NormalizedCellType t=(INTERP_KERNEL::NormalizedCellType)*conn;
+      switch(t)
+        {
+          case INTERP_KERNEL::NORM_QUAD4:
+            {
+              fillInCompact3DMode(3,4,conn+1,coo,tmp);
+              *pt=INTERP_KERNEL::quadSkew(tmp);
+              break;
+            }
+        default:
+          throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getSkewField : A cell with not manged type (NORM_QUAD4) has been detected !");
+        }
+      conn+=connI[i+1]-connI[i];
+    }
+  ret->setName("Skew");
+  ret->incrRef();
+  return ret;
+}
+
 /*!
  * This method aggregate the bbox of each cell and put it into bbox parameter.
  * @param bbox out parameter of size 2*spacedim*nbOfcells.
@@ -2891,3 +3108,25 @@ void MEDCouplingUMesh::tryToCorrectPolyhedronOrientation(int *begin, int *end, c
         }
     }
 }
+
+/*!
+ * This method put in zip format into parameter 'zipFrmt' in full interlace mode.
+ * This format is often asked by INTERP_KERNEL algorithms to avoid many indirections into coordinates array.
+ */
+void MEDCouplingUMesh::fillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt) throw(INTERP_KERNEL::Exception)
+{
+  double *w=zipFrmt;
+  if(spaceDim==3)
+    for(int i=0;i<nbOfNodesInCell;i++)
+      w=std::copy(coo+3*conn[i],coo+3*conn[i]+3,w);
+  else if(spaceDim==2)
+    {
+      for(int i=0;i<nbOfNodesInCell;i++)
+        {
+          w=std::copy(coo+2*conn[i],coo+2*conn[i]+2,w);
+          *w++=0.;
+        }
+    }
+  else
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::fillInCompact3DMode : Invalid spaceDim specified : must be 2 or 3 !");
+}
diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 8c4f99882..eb989c9a6 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -113,6 +113,11 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT void arePolyhedronsNotCorrectlyOriented(std::vector<int>& cells) const throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT void getFastAveragePlaneOfThis(double *vec, double *pos) const throw(INTERP_KERNEL::Exception);
+    //Mesh quality
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getEdgeRatioField() const throw(INTERP_KERNEL::Exception);
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getAspectRatioField() const throw(INTERP_KERNEL::Exception);
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getWarpField() const throw(INTERP_KERNEL::Exception);
+    MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getSkewField() const throw(INTERP_KERNEL::Exception);
     //utilities for MED File RW
     MEDCOUPLING_EXPORT bool checkConsecutiveCellTypes() const;
     MEDCOUPLING_EXPORT bool checkConsecutiveCellTypesAndOrder(const INTERP_KERNEL::NormalizedCellType *orderBg, const INTERP_KERNEL::NormalizedCellType *orderEnd) const;
@@ -148,6 +153,7 @@ namespace ParaMEDMEM
     template<int SPACEDIM>
     void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
                                      double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const;
+    static void fillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt) throw(INTERP_KERNEL::Exception);
     static void appendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector<int>& ret);
   private:
     //! this iterator stores current position in _nodal_connec array.
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index 9bbdb8519..8382a9230 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -124,6 +124,7 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( testSortPerTuple1 );
     CPPUNIT_TEST( testIsEqualWithoutConsideringStr1 );
     CPPUNIT_TEST( testGetNodeIdsOfCell1 );
+    CPPUNIT_TEST( testGetEdgeRatioField1 );
     //MEDCouplingBasicsTestInterp.cxx
     CPPUNIT_TEST( test2DInterpP0P0_1 );
     CPPUNIT_TEST( test2DInterpP0P0PL_1 );
@@ -276,6 +277,7 @@ namespace ParaMEDMEM
     void testSortPerTuple1();
     void testIsEqualWithoutConsideringStr1();
     void testGetNodeIdsOfCell1();
+    void testGetEdgeRatioField1();
     //MEDCouplingBasicsTestInterp.cxx
     void test2DInterpP0P0_1();
     void test2DInterpP0P0PL_1();
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx
index a06946267..f636313cf 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx
@@ -1977,3 +1977,28 @@ void MEDCouplingBasicsTest::testGetNodeIdsOfCell1()
   CPPUNIT_ASSERT_DOUBLES_EQUAL(0.2,coords[1],1e-13);
   mesh1->decrRef();
 }
+
+void MEDCouplingBasicsTest::testGetEdgeRatioField1()
+{
+  MEDCouplingUMesh *m1=build2DTargetMesh_1();
+  MEDCouplingFieldDouble *f1=m1->getEdgeRatioField();
+  CPPUNIT_ASSERT_EQUAL(m1->getNumberOfCells(),f1->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(5,f1->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(1,f1->getNumberOfComponents());
+  const double expected1[5]={1.,1.4142135623730951, 1.4142135623730951,1.,1.};
+  for(int i=0;i<5;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected1[i],f1->getIJ(i,0),1e-14);
+  f1->decrRef();
+  m1->decrRef();
+  //
+  m1=build3DSurfTargetMesh_1();
+  f1=m1->getEdgeRatioField();
+  CPPUNIT_ASSERT_EQUAL(m1->getNumberOfCells(),f1->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(5,f1->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(1,f1->getNumberOfComponents());
+  const double expected2[5]={1.4142135623730951, 1.7320508075688772, 1.7320508075688772, 1.4142135623730951, 1.4142135623730951};
+  for(int i=0;i<5;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected2[i],f1->getIJ(i,0),1e-14);
+  f1->decrRef();
+  m1->decrRef();
+}
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 6120fa854..53cd9f4df 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -3377,6 +3377,28 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         li2=mesh1.getNodalConnectivityIndex().getValuesAsTuple()
         self.assertEqual(6,len(li2))
         pass
+
+    def testGetEdgeRatioField1(self):
+        m1=MEDCouplingDataForTest.build2DTargetMesh_1();
+        f1=m1.getEdgeRatioField();
+        self.assertEqual(m1.getNumberOfCells(),f1.getNumberOfTuples());
+        self.assertEqual(5,f1.getNumberOfTuples());
+        self.assertEqual(1,f1.getNumberOfComponents());
+        expected1=[1.,1.4142135623730951, 1.4142135623730951,1.,1.]
+        for i in xrange(5):
+            self.assertAlmostEqual(expected1[i],f1.getIJ(i,0),14);
+            pass
+        #
+        m1=MEDCouplingDataForTest.build3DSurfTargetMesh_1();
+        f1=m1.getEdgeRatioField();
+        self.assertEqual(m1.getNumberOfCells(),f1.getNumberOfTuples());
+        self.assertEqual(5,f1.getNumberOfTuples());
+        self.assertEqual(1,f1.getNumberOfComponents());
+        expected2=[1.4142135623730951, 1.7320508075688772, 1.7320508075688772, 1.4142135623730951, 1.4142135623730951]
+        for i in xrange(5):
+            self.assertAlmostEqual(expected2[i],f1.getIJ(i,0),14);
+            pass
+        pass
     
     def setUp(self):
         pass
diff --git a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
index 09aef4ff4..48a7115d9 100644
--- a/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
+++ b/src/MEDCoupling_Swig/libMEDCoupling_Swig.i
@@ -141,6 +141,10 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::convertCellArrayPerGeoType;
 %newobject ParaMEDMEM::MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildDirectionVectorField;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getEdgeRatioField;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getAspectRatioField;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getWarpField;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getSkewField;
 %newobject ParaMEDMEM::MEDCouplingExtrudedMesh::New;
 %newobject ParaMEDMEM::MEDCouplingExtrudedMesh::build3DUnstructuredMesh;
 %newobject ParaMEDMEM::MEDCouplingCMesh::New;
@@ -479,6 +483,10 @@ namespace ParaMEDMEM
     bool isPresenceOfQuadratic() const;
     MEDCouplingFieldDouble *buildDirectionVectorField() const;
     void convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception);
+    MEDCouplingFieldDouble *getEdgeRatioField() const throw(INTERP_KERNEL::Exception);
+    MEDCouplingFieldDouble *getAspectRatioField() const throw(INTERP_KERNEL::Exception);
+    MEDCouplingFieldDouble *getWarpField() const throw(INTERP_KERNEL::Exception);
+    MEDCouplingFieldDouble *getSkewField() const throw(INTERP_KERNEL::Exception);
     %extend {
       std::string __str__() const
       {