From: Anthony Geay Date: Wed, 29 Apr 2015 06:58:30 +0000 (+0200) Subject: Implementation of MEDCouplingPointSet.computeDiameterField (EDF10718) X-Git-Tag: V7_6_0rc1~6 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=refs%2Fheads%2Fagy%2FComputeDiameter;p=tools%2Fmedcoupling.git Implementation of MEDCouplingPointSet.computeDiameterField (EDF10718) --- diff --git a/src/INTERP_KERNEL/CMakeLists.txt b/src/INTERP_KERNEL/CMakeLists.txt index 3523c9b85..1fe07853b 100644 --- a/src/INTERP_KERNEL/CMakeLists.txt +++ b/src/INTERP_KERNEL/CMakeLists.txt @@ -26,6 +26,7 @@ SET(interpkernel_SOURCES TranslationRotationMatrix.cxx TetraAffineTransform.cxx CellModel.cxx + DiameterCalculator.cxx UnitTetraIntersectionBary.cxx InterpolationOptions.cxx BoxSplittingOptions.cxx diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx index 959162b0a..7c01d990c 100644 --- a/src/INTERP_KERNEL/CellModel.cxx +++ b/src/INTERP_KERNEL/CellModel.cxx @@ -21,6 +21,7 @@ #include "CellModel.hxx" #include "InterpKernelException.hxx" +#include "DiameterCalculator.hxx" #include #include @@ -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 !"); + } + } } diff --git a/src/INTERP_KERNEL/CellModel.hxx b/src/INTERP_KERNEL/CellModel.hxx index 6b709012e..d1b5e1d95 100644 --- a/src/INTERP_KERNEL/CellModel.hxx +++ b/src/INTERP_KERNEL/CellModel.hxx @@ -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 index 000000000..a27395996 --- /dev/null +++ b/src/INTERP_KERNEL/DiameterCalculator.cxx @@ -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 +#include +#include + +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 +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 +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 +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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorTRI3S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorTRI3S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorQUAD4S2::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorQUAD4S3::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorTETRA4::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorHEXA8::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorPENTA6::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(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(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(bgId,endId,indPtr,connPtr,coordsPtr,resPtr); +} + +void DiameterCalulatorPYRA5::computeFor1SGTUMeshFrmt(int nbOfCells, const int *connPtr, const double *coordsPtr, double *resPtr) const +{ + ComputeFor1SGTUMeshFrmt(nbOfCells,connPtr,coordsPtr,resPtr); +} diff --git a/src/INTERP_KERNEL/DiameterCalculator.hxx b/src/INTERP_KERNEL/DiameterCalculator.hxx new file mode 100644 index 000000000..865e86682 --- /dev/null +++ b/src/INTERP_KERNEL/DiameterCalculator.hxx @@ -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 diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx index 343076dfd..eb0db7da8 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.cxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.cxx @@ -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 ret(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME)); + int nbCells(getNumberOfCells()); + MEDCouplingAutoRefCountObjectPtr arr(DataArrayDouble::New()); + arr->alloc(nbCells,1); + INTERP_KERNEL::AutoCppPtr 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 MEDCoupling1DGTUMesh::BuildAPolygonFromParts(const std::vector< std::vector >& parts) { std::vector ret; diff --git a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx index 3ccb85cb5..9928bc48f 100644 --- a/src/MEDCoupling/MEDCoupling1GTUMesh.hxx +++ b/src/MEDCoupling/MEDCoupling1GTUMesh.hxx @@ -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); diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index 08281ec96..21c830d2f 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -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; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index f213a9d42..b8421f251 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -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 ret(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME)); + ret->setMesh(this); + std::set types; + ComputeAllTypesInternal(types,_nodal_connec,_nodal_connec_index); + int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells()); + MEDCouplingAutoRefCountObjectPtr arr(DataArrayDouble::New()); + arr->alloc(nbCells,1); + for(std::set::const_iterator it=types.begin();it!=types.end();it++) + { + INTERP_KERNEL::AutoCppPtr dc(INTERP_KERNEL::CellModel::GetCellModel(*it).buildInstanceOfDiameterCalulator(spaceDim)); + MEDCouplingAutoRefCountObjectPtr 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& 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) { diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index a3972a295..fff81f7c5 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -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 getDistributionOfTypes() const; MEDCOUPLING_EXPORT DataArrayInt *checkTypeConsistencyAndContig(const std::vector& code, const std::vector& 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& types, const DataArrayInt *nodalConnec, const DataArrayInt *nodalConnecIndex); public: MEDCOUPLING_EXPORT static DataArrayInt *ComputeRangesFromTypeDistribution(const std::vector& code); MEDCOUPLING_EXPORT static const int N_MEDMEM_ORDER=24; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index c6bc0d7d0..bcfa2c863 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -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__': diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index 6926805ae..22cd18921 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -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)