From 98db8972a3182e1ef4d849e9745d54d5682696bd Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 18 Oct 2010 07:43:53 +0000 Subject: [PATCH] *** empty log message *** --- src/MEDCoupling/MEDCouplingUMesh.cxx | 239 ++++++++++++++++++ src/MEDCoupling/MEDCouplingUMesh.hxx | 6 + .../Test/MEDCouplingBasicsTest.hxx | 2 + .../Test/MEDCouplingBasicsTest2.cxx | 25 ++ src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 22 ++ src/MEDCoupling_Swig/libMEDCoupling_Swig.i | 8 + 6 files changed, 302 insertions(+) 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 @@ -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 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;isetName("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 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;isetName("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 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;isetName("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 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;isetName("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& 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 void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, double eps, std::vector& elts, std::vector& 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& 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 { -- 2.39.2