From 43f01d39f46248e7c5ff1c6b3665f11ab5a7c5b4 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Tue, 20 Dec 2022 16:13:16 +0100 Subject: [PATCH] Implementation of MEDCoupling1SGTUMesh.computeTriangleHeight and DataArrayDouble.minPerTuple --- src/INTERP_KERNEL/VolSurfUser.cxx | 13 ++++- src/INTERP_KERNEL/VolSurfUser.hxx | 3 ++ src/MEDCoupling/MEDCoupling1GTUMesh.cxx | 47 +++++++++++++++++++ src/MEDCoupling/MEDCoupling1GTUMesh.hxx | 1 + src/MEDCoupling/MEDCouplingMemArray.cxx | 41 +++++++++++----- src/MEDCoupling/MEDCouplingMemArray.hxx | 4 ++ .../MEDCouplingBasicsTest7.py | 13 +++++ src/MEDCoupling_Swig/MEDCouplingCommon.i | 7 +++ src/MEDCoupling_Swig/MEDCouplingMemArray.i | 2 + 9 files changed, 119 insertions(+), 12 deletions(-) diff --git a/src/INTERP_KERNEL/VolSurfUser.cxx b/src/INTERP_KERNEL/VolSurfUser.cxx index 924d2bef5..fb9bb90f0 100644 --- a/src/INTERP_KERNEL/VolSurfUser.cxx +++ b/src/INTERP_KERNEL/VolSurfUser.cxx @@ -21,6 +21,7 @@ #include "VolSurfUser.hxx" #include "InterpKernelAutoPtr.hxx" #include "InterpolationUtils.hxx" +#include "VectorUtils.hxx" #include #include @@ -276,5 +277,15 @@ namespace INTERP_KERNEL matrix[11]=-p0[0]*matrix[8]-p0[1]*matrix[9]-p0[2]*matrix[10]; return true; } -} + template + void ComputeTriangleHeight(const double *PA, const double *PB, const double *PC, double *res) + { + double AB = getDistanceBtw2Pts(PA,PB); + double BC = getDistanceBtw2Pts(PB,PC); + double CA = getDistanceBtw2Pts(PC,PA); + double perim( (AB+BC+CA)*0.5 ); + double num( 2*sqrt(perim*(perim-AB)*(perim-BC)*(perim-CA)) ); + res[0] = num/AB; res[1] = num/BC; res[2] = num/CA; + } +} diff --git a/src/INTERP_KERNEL/VolSurfUser.hxx b/src/INTERP_KERNEL/VolSurfUser.hxx index 79e346f08..4682ed7ed 100644 --- a/src/INTERP_KERNEL/VolSurfUser.hxx +++ b/src/INTERP_KERNEL/VolSurfUser.hxx @@ -49,6 +49,9 @@ namespace INTERP_KERNEL double INTERPKERNEL_EXPORT DistanceFromPtToPolygonInSpaceDim3(const double *pt, const mcIdType *connOfPolygonBg, const mcIdType *connOfPolygonEnd, const double *coords); bool ComputeRotTranslationMatrixToPut3PointsOnOXY(const double *pt0Tri3, const double *pt1Tri3, const double *pt2Tri3, double *matrix); + + template + void ComputeTriangleHeight(const double *PA, const double *PB, const double *PC, double *res); } #endif diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx index d0a5917c3..8b6567583 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx @@ -27,6 +27,7 @@ #include "DiameterCalculator.hxx" #include "OrientationInverter.hxx" #include "InterpKernelAutoPtr.hxx" +#include "VolSurfUser.cxx" using namespace MEDCoupling; @@ -1743,6 +1744,52 @@ MEDCoupling1SGTUMesh *MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4() const return ret.retn(); } +/*! + * This method for each cell in \a this the triangle height for each edge in a newly allocated/created array instance. + * + * \return DataArrayDouble * - a newly allocated instance with this->getNumberOfCells() tuples and 3 components storing for each cell in \a this the corresponding height. + * \throw If \a this is not a mesh containing only NORM_TRI3 cells. + * \throw If \a this is not properly allocated. + * \throw If spaceDimension is not in 2 or 3. + */ +MCAuto MEDCoupling1SGTUMesh::computeTriangleHeight() const +{ + checkConsistencyLight(); + const INTERP_KERNEL::CellModel& cm(getCellModel()); + if(cm.getEnum()!=INTERP_KERNEL::NORM_TRI3) + THROW_IK_EXCEPTION("MEDCoupling1SGTUMesh::computeTriangleHeight : this method can be applied only on TRI3 mesh !"); + MCAuto ret(DataArrayDouble::New()); + mcIdType nbTri3( getNumberOfCells() ); + const double *coordPtr( this->getCoords()->begin() ); + const mcIdType *inConnPtr(getNodalConnectivity()->begin()); + ret->alloc(nbTri3,3); + double *retPtr( ret->getPointer() ); + switch( this->getSpaceDimension()) + { + case 2: + { + constexpr unsigned SPACEDIM = 2; + for(mcIdType iCell = 0 ; iCell < nbTri3 ; ++iCell) + { + INTERP_KERNEL::ComputeTriangleHeight(coordPtr + SPACEDIM*inConnPtr[3*iCell+0], coordPtr + SPACEDIM*inConnPtr[3*iCell+1], coordPtr + SPACEDIM*inConnPtr[3*iCell+2],retPtr+3*iCell); + } + break; + } + case 3: + { + constexpr unsigned SPACEDIM = 3; + for(mcIdType iCell = 0 ; iCell < nbTri3 ; ++iCell) + { + INTERP_KERNEL::ComputeTriangleHeight(coordPtr + SPACEDIM*inConnPtr[3*iCell+0], coordPtr + SPACEDIM*inConnPtr[3*iCell+1], coordPtr + SPACEDIM*inConnPtr[3*iCell+2],retPtr+3*iCell); + } + break; + } + default: + THROW_IK_EXCEPTION("MEDCoupling1SGTUMesh::computeTriangleHeight : only spacedim in [2,3] supported !"); + } + return ret; +} + /*! * This method starts from an unstructured mesh that hides in reality a cartesian mesh. * If it is not the case, an exception will be thrown. diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx index c5ee4fb32..fdfc7a745 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx @@ -158,6 +158,7 @@ namespace MEDCoupling MEDCOUPLING_EXPORT MEDCoupling1GTUMesh *computeDualMesh() const; MEDCOUPLING_EXPORT DataArrayIdType *sortHexa8EachOther(); MEDCOUPLING_EXPORT MEDCoupling1SGTUMesh *explodeEachHexa8To6Quad4() const; + MEDCOUPLING_EXPORT MCAuto computeTriangleHeight() const; MEDCOUPLING_EXPORT MEDCouplingCMesh *structurizeMe(DataArrayIdType *& cellPerm, DataArrayIdType *& nodePerm, double eps=1e-12) const; public://serialization MEDCOUPLING_EXPORT void getTinySerializationInformation(std::vector& tinyInfoD, std::vector& tinyInfo, std::vector& littleStrings) const; diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 60d040806..5006c10a8 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -2337,16 +2337,7 @@ DataArrayDouble *DataArrayDouble::magnitude() const return ret; } -/*! - * Computes the maximal value within every tuple of \a this array. - * \return DataArrayDouble * - the new instance of DataArrayDouble containing the - * same number of tuples as \a this array and one component. - * The caller is to delete this result array using decrRef() as it is no more - * needed. - * \throw If \a this is not allocated. - * \sa DataArrayDouble::maxPerTupleWithCompoId - */ -DataArrayDouble *DataArrayDouble::maxPerTuple() const +DataArrayDouble *DataArrayDouble::operatePerTuple(std::function func) const { checkAllocated(); std::size_t nbOfComp(getNumberOfComponents()); @@ -2356,10 +2347,38 @@ DataArrayDouble *DataArrayDouble::maxPerTuple() const const double *src=getConstPointer(); double *dest=ret->getPointer(); for(mcIdType i=0;ioperatePerTuple([](const double *bg, const double *endd) { return *std::max_element(bg,endd); }); +} + +/*! + * Computes the minimal value within every tuple of \a this array. + * \return DataArrayDouble * - the new instance of DataArrayDouble containing the + * same number of tuples as \a this array and one component. + * The caller is to delete this result array using decrRef() as it is no more + * needed. + * \throw If \a this is not allocated. + * \sa DataArrayDouble::maxPerTuple + */ +DataArrayDouble *DataArrayDouble::minPerTuple() const +{ + return this->operatePerTuple([](const double *bg, const double *endd) { return *std::min_element(bg,endd); }); +} + /*! * Computes the maximal value within every tuple of \a this array and it returns the first component * id for each tuple that corresponds to the maximal value within the tuple. diff --git a/src/MEDCoupling/MEDCouplingMemArray.hxx b/src/MEDCoupling/MEDCouplingMemArray.hxx index cc846eeae..ea18707e3 100755 --- a/src/MEDCoupling/MEDCouplingMemArray.hxx +++ b/src/MEDCoupling/MEDCouplingMemArray.hxx @@ -33,6 +33,7 @@ #include #include #include +#include namespace MEDCoupling { @@ -497,6 +498,7 @@ namespace MEDCoupling DataArrayDouble *trace() const; DataArrayDouble *deviator() const; DataArrayDouble *magnitude() const; + DataArrayDouble *minPerTuple() const; DataArrayDouble *maxPerTuple() const; DataArrayDouble *maxPerTupleWithCompoId(DataArrayIdType* &compoIdOfMaxPerTuple) const; DataArrayDouble *buildEuclidianDistanceDenseMatrix() const; @@ -544,6 +546,8 @@ namespace MEDCoupling template static void FindTupleIdsNearTuplesAlg(const BBTreePts& myTree, const double *pos, mcIdType nbOfTuples, double eps, DataArrayIdType *c, DataArrayIdType *cI); + private: + DataArrayDouble *operatePerTuple(std::function func) const; private: ~DataArrayDouble() { } DataArrayDouble() { } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index 36d144b4f..a9a0c6006 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1119,6 +1119,19 @@ class MEDCouplingBasicsTest7(unittest.TestCase): delta_Z.abs() self.assertTrue(delta_Z.findIdsNotInRange(-1e-5,+1e-5).empty()) + def testComputeTriangleHeight0(self): + arr = DataArrayDouble([0,1]) + m = MEDCouplingCMesh() ; m.setCoords(arr,arr) + m = m.buildUnstructured() ; m.simplexize(0) ; m = MEDCoupling1SGTUMesh(m) + res = m.computeTriangleHeight() + expected = DataArrayDouble([(1.0, 1.0, sqrt(2)/2.0), (sqrt(2)/2.0, 1.0, 1.0)]) + self.assertTrue( res.isEqual(expected,1e-12) ) + m.changeSpaceDimension(3,100) + res2 = m.computeTriangleHeight() + self.assertTrue( res2.isEqual(expected,1e-12) ) + expected2 = DataArrayDouble([sqrt(2)/2.0, sqrt(2)/2.0]) + self.assertTrue( res2.minPerTuple().isEqual(expected2,1e-12) ) + pass if __name__ == '__main__': diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 7acf25f32..1793504a5 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -413,6 +413,7 @@ typedef long mcPyPtrType; %newobject MEDCoupling::MEDCoupling1SGTUMesh::buildSetInstanceFromThis; %newobject MEDCoupling::MEDCoupling1SGTUMesh::computeDualMesh; %newobject MEDCoupling::MEDCoupling1SGTUMesh::explodeEachHexa8To6Quad4; +%newobject MEDCoupling::MEDCoupling1SGTUMesh::computeTriangleHeight; %newobject MEDCoupling::MEDCoupling1SGTUMesh::sortHexa8EachOther; %newobject MEDCoupling::MEDCoupling1SGTUMesh::Merge1SGTUMeshes; %newobject MEDCoupling::MEDCoupling1SGTUMesh::Merge1SGTUMeshesOnSameCoords; @@ -3204,6 +3205,12 @@ namespace MEDCoupling return ret; } + DataArrayDouble *MEDCoupling1SGTUMesh::computeTriangleHeight() const + { + MCAuto ret = self->computeTriangleHeight(); + return ret.retn(); + } + static MEDCoupling1SGTUMesh *Merge1SGTUMeshes(PyObject *li) { std::vector tmp; diff --git a/src/MEDCoupling_Swig/MEDCouplingMemArray.i b/src/MEDCoupling_Swig/MEDCouplingMemArray.i index bacf22e55..f9ba14216 100644 --- a/src/MEDCoupling_Swig/MEDCouplingMemArray.i +++ b/src/MEDCoupling_Swig/MEDCouplingMemArray.i @@ -288,6 +288,7 @@ %newobject MEDCoupling::DataArrayDouble::deviator; %newobject MEDCoupling::DataArrayDouble::magnitude; %newobject MEDCoupling::DataArrayDouble::maxPerTuple; +%newobject MEDCoupling::DataArrayDouble::minPerTuple; %newobject MEDCoupling::DataArrayDouble::sumPerTuple; %newobject MEDCoupling::DataArrayDouble::computeBBoxPerTuple; %newobject MEDCoupling::DataArrayDouble::buildEuclidianDistanceDenseMatrix; @@ -1048,6 +1049,7 @@ typedef DataArrayInt64 DataArrayIdType; DataArrayDouble *deviator() const; DataArrayDouble *magnitude() const; DataArrayDouble *maxPerTuple() const; + DataArrayDouble *minPerTuple() const; DataArrayDouble *sumPerTuple() const; DataArrayDouble *buildEuclidianDistanceDenseMatrix() const; DataArrayDouble *buildEuclidianDistanceDenseMatrixWith(const DataArrayDouble *other) const; -- 2.39.2