From: ageay Date: Fri, 28 Aug 2009 08:30:59 +0000 (+0000) Subject: *** empty log message *** X-Git-Tag: V5_1_main_FINAL~353 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=8909ba90cc355552bbb94770942f6d833675b987;p=tools%2Fmedcoupling.git *** empty log message *** --- diff --git a/src/INTERP_KERNEL/IntegralUniformIntersector.hxx b/src/INTERP_KERNEL/IntegralUniformIntersector.hxx new file mode 100644 index 000000000..8297323e6 --- /dev/null +++ b/src/INTERP_KERNEL/IntegralUniformIntersector.hxx @@ -0,0 +1,70 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTEGRALUNIFORMINTERSECTOR_HXX__ +#define __INTEGRALUNIFORMINTERSECTOR_HXX__ + +#include "TargetIntersector.hxx" + +#include + +namespace INTERP_KERNEL +{ + template + class IntegralUniformIntersector : public TargetIntersector + { + public: + typedef typename MyMeshType::MyConnType ConnType; + public: + IntegralUniformIntersector(const MyMeshType& mesh, bool isAbs); + double performNormalization(double val) const { if(_is_abs) return fabs(val); else return val; } + void setFromTo(bool val) { _from_to=val; } + void putValueIn(ConnType i, double val, MyMatrix& res) const; + protected: + const MyMeshType& _mesh; + //! if false means fromIntegralUniform if true means toIntegralUniform + bool _from_to; + bool _is_abs; + }; + + template + class IntegralUniformIntersectorP0 : public IntegralUniformIntersector + { + public: + typedef typename MyMeshType::MyConnType ConnType; + public: + IntegralUniformIntersectorP0(const MyMeshType& mesh, bool isAbs); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + void intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res); + }; + + template + class IntegralUniformIntersectorP1 : public IntegralUniformIntersector + { + public: + typedef typename MyMeshType::MyConnType ConnType; + public: + IntegralUniformIntersectorP1(const MyMeshType& mesh, bool isAbs); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + void intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/IntegralUniformIntersector.txx b/src/INTERP_KERNEL/IntegralUniformIntersector.txx new file mode 100644 index 000000000..0acc2a18c --- /dev/null +++ b/src/INTERP_KERNEL/IntegralUniformIntersector.txx @@ -0,0 +1,155 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTEGRALUNIFORMINTERSECTOR_TXX__ +#define __INTEGRALUNIFORMINTERSECTOR_TXX__ + +#include "IntegralUniformIntersector.hxx" +#include "VolSurfUser.txx" + +namespace INTERP_KERNEL +{ + template + IntegralUniformIntersector::IntegralUniformIntersector(const MyMeshType& mesh, bool isAbs):_mesh(mesh),_from_to(false),_is_abs(isAbs) + { + } + + template + void IntegralUniformIntersector::putValueIn(ConnType iInCMode, double val1, MyMatrix& res) const + { + static const NumberingPolicy numPol=MyMeshType::My_numPol; + double val=performNormalization(val1); + if(_from_to) + { + typename MyMatrix::value_type& resRow=res[0]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(iInCMode)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(iInCMode),val)); + else + { + double val2=(*iterRes).second+val; + resRow.erase(OTT::indFC(iInCMode)); + resRow.insert(std::make_pair(OTT::indFC(iInCMode),val2)); + } + } + else + { + typename MyMatrix::value_type& resRow=res[iInCMode]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(0)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(0),val)); + else + { + double val2=(*iterRes).second+val; + resRow.erase(OTT::indFC(0)); + resRow.insert(std::make_pair(OTT::indFC(0),val2)); + } + } + } + + template + IntegralUniformIntersectorP0::IntegralUniformIntersectorP0(const MyMeshType& mesh, bool isAbs):IntegralUniformIntersector(mesh,isAbs) + { + } + + template + int IntegralUniformIntersectorP0::getNumberOfRowsOfResMatrix() const + { + if(IntegralUniformIntersector::_from_to) + return 1; + else + return IntegralUniformIntersector::_mesh.getNumberOfElements(); + } + + template + int IntegralUniformIntersectorP0::getNumberOfColsOfResMatrix() const + { + if(IntegralUniformIntersector::_from_to) + return IntegralUniformIntersector::_mesh.getNumberOfElements(); + else + return 1; + } + + template + void IntegralUniformIntersectorP0::intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res) + { + static const NumberingPolicy numPol=MyMeshType::My_numPol; + res.resize(getNumberOfRowsOfResMatrix()); + unsigned long nbelem=IntegralUniformIntersector::_mesh.getNumberOfElements(); + const ConnType *connIndx=IntegralUniformIntersector::_mesh.getConnectivityIndexPtr(); + const ConnType *conn=IntegralUniformIntersector::_mesh.getConnectivityPtr(); + const double *coords=IntegralUniformIntersector::_mesh.getCoordinatesPtr(); + for(unsigned long i=0;i::_mesh.getTypeOfElement(OTT::indFC(i)); + double val=computeVolSurfOfCell(t,conn+OTT::ind2C(connIndx[i]),connIndx[i+1]-connIndx[i],coords); + IntegralUniformIntersector::putValueIn(i,val,res); + } + } + + template + IntegralUniformIntersectorP1::IntegralUniformIntersectorP1(const MyMeshType& mesh, bool isAbs):IntegralUniformIntersector(mesh,isAbs) + { + } + + template + int IntegralUniformIntersectorP1::getNumberOfRowsOfResMatrix() const + { + if(IntegralUniformIntersector::_from_to) + return 1; + else + return IntegralUniformIntersector::_mesh.getNumberOfNodes(); + } + + template + int IntegralUniformIntersectorP1::getNumberOfColsOfResMatrix() const + { + if(IntegralUniformIntersector::_from_to) + return IntegralUniformIntersector::_mesh.getNumberOfNodes(); + else + return 1; + } + + template + void IntegralUniformIntersectorP1::intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res) + { + static const NumberingPolicy numPol=MyMeshType::My_numPol; + res.resize(getNumberOfRowsOfResMatrix()); + unsigned long nbelem=IntegralUniformIntersector::_mesh.getNumberOfElements(); + const ConnType *connIndx=IntegralUniformIntersector::_mesh.getConnectivityIndexPtr(); + const ConnType *conn=IntegralUniformIntersector::_mesh.getConnectivityPtr(); + const double *coords=IntegralUniformIntersector::_mesh.getCoordinatesPtr(); + for(unsigned long i=0;i::_mesh.getTypeOfElement(OTT::indFC(i)); + int lgth=connIndx[i+1]-connIndx[i]; + const ConnType *locConn=conn+OTT::ind2C(connIndx[i]); + double val=computeVolSurfOfCell(t,locConn,lgth,coords); + if(t==NORM_TRI3) + val/=3.; + else if(t==NORM_TETRA4) + val/=4.; + else + throw INTERP_KERNEL::Exception("Invalid cell type detected : must be TRI3 or TETRA4 ! "); + for(int j=0;j::putValueIn(OTT::coo2C(locConn[j]),val,res); + } + } +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx index 57edda88e..fe34d1b02 100644 --- a/src/INTERP_KERNEL/Interpolation.hxx +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -38,8 +38,15 @@ namespace INTERP_KERNEL Interpolation(const InterpolationOptions& io) :InterpolationOptions(io){} //interpolation of two triangular meshes. template - int interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result) - { return asLeaf().interpolateMeshes(mesh1,mesh2,result); } + int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result) + { return asLeaf().interpolateMeshes(meshS,meshT,result); } + template + int fromIntegralUniform(const MyMeshType& meshT, MatrixType& result, const char *method) { return fromToIntegralUniform(true,meshT,result,method); } + template + int toIntegralUniform(const MyMeshType& meshS, MatrixType& result, const char *method) { return fromToIntegralUniform(false,meshS,result,method); } + protected: + template + int fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method); protected: TrueMainInterpolator& asLeaf() { return static_cast(*this); } }; diff --git a/src/INTERP_KERNEL/Interpolation.txx b/src/INTERP_KERNEL/Interpolation.txx new file mode 100644 index 000000000..483d8d1e8 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation.txx @@ -0,0 +1,58 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __INTERPOLATION_TXX__ +#define __INTERPOLATION_TXX__ + +#include "Interpolation.hxx" +#include "IntegralUniformIntersector.hxx" +#include "IntegralUniformIntersector.txx" + +namespace INTERP_KERNEL +{ + template + template + int Interpolation::fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method) + { + typedef typename MyMeshType::MyConnType ConnType; + std::string methodCPP(method); + int ret=-1; + if(methodCPP=="P0") + { + IntegralUniformIntersectorP0 intersector(mesh,InterpolationOptions::getMeasureAbsStatus()); + intersector.setFromTo(fromTo); + std::vector tmp; + intersector.intersectCells(0,tmp,result); + ret=intersector.getNumberOfColsOfResMatrix(); + } + else if(methodCPP=="P1") + { + IntegralUniformIntersectorP1 intersector(mesh,InterpolationOptions::getMeasureAbsStatus()); + intersector.setFromTo(fromTo); + std::vector tmp; + intersector.intersectCells(0,tmp,result); + ret=intersector.getNumberOfColsOfResMatrix(); + } + else + throw INTERP_KERNEL::Exception("Invalid method specified in fromIntegralUniform : must be in { \"P0\", \"P1\"}"); + return ret; + } +} + +#endif + diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx index d24cbdc8a..35ee63b34 100644 --- a/src/INTERP_KERNEL/Interpolation3D.txx +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -20,6 +20,7 @@ #define __INTERPOLATION3D_TXX__ #include "Interpolation3D.hxx" +#include "Interpolation.txx" #include "MeshElement.txx" #include "TransformedTriangle.hxx" #include "PolyhedronIntersector.txx" @@ -315,7 +316,6 @@ namespace INTERP_KERNEL return ret; } - } #endif diff --git a/src/INTERP_KERNEL/InterpolationOptions.cxx b/src/INTERP_KERNEL/InterpolationOptions.cxx index 6bbb520b6..c4514b160 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.cxx +++ b/src/INTERP_KERNEL/InterpolationOptions.cxx @@ -40,6 +40,8 @@ const char INTERP_KERNEL::InterpolationOptions::DO_ROTATE_STR[]="DoRotate"; const char INTERP_KERNEL::InterpolationOptions::ORIENTATION_STR[]="Orientation"; +const char INTERP_KERNEL::InterpolationOptions::MEASURE_ABS_STR[]="MeasureAbs"; + const char INTERP_KERNEL::InterpolationOptions::INTERSEC_TYPE_STR[]="IntersectionType"; const char INTERP_KERNEL::InterpolationOptions::SPLITTING_POLICY_STR[]="SplittingPolicy"; @@ -106,6 +108,11 @@ bool INTERP_KERNEL::InterpolationOptions::setOptionInt(const std::string& key, i setOrientation(value); return true; } + else if(key==MEASURE_ABS_STR) + { + bool valBool=(value!=0); + setMeasureAbsStatus(valBool); + } else return false; } diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx index bfc843508..89430d1f0 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.hxx +++ b/src/INTERP_KERNEL/InterpolationOptions.hxx @@ -46,6 +46,7 @@ namespace INTERP_KERNEL double _bounding_box_adjustment_abs ; double _max_distance_for_3Dsurf_intersect; int _orientation ; + bool _measure_abs; SplittingPolicy _splitting_policy ; public: @@ -76,6 +77,9 @@ namespace INTERP_KERNEL int getOrientation() const { return _orientation; } void setOrientation(int o) { _orientation=o; } + + bool getMeasureAbsStatus() const { return _measure_abs; } + void setMeasureAbsStatus(bool newStatus) { _measure_abs=newStatus; } SplittingPolicy getSplittingPolicy() const { return _splitting_policy; } void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; } @@ -90,6 +94,7 @@ namespace INTERP_KERNEL _bounding_box_adjustment_abs=0.; _max_distance_for_3Dsurf_intersect=DFT_MAX_DIST_3DSURF_INTERSECT; _orientation=0; + _measure_abs=true; _splitting_policy=GENERAL_48; } bool setOptionDouble(const std::string& key, double value); @@ -108,6 +113,7 @@ namespace INTERP_KERNEL static const char PRINT_LEV_STR[]; static const char DO_ROTATE_STR[]; static const char ORIENTATION_STR[]; + static const char MEASURE_ABS_STR[]; static const char INTERSEC_TYPE_STR[]; static const char SPLITTING_POLICY_STR[]; static const char TRIANGULATION_INTERSECT2D_STR[]; diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx index f4321c95e..51bd62773 100755 --- a/src/INTERP_KERNEL/InterpolationPlanar.hxx +++ b/src/INTERP_KERNEL/InterpolationPlanar.hxx @@ -31,8 +31,6 @@ namespace INTERP_KERNEL { private: double _dim_caracteristic; - static const double DEFAULT_PRECISION; - public: InterpolationPlanar(); InterpolationPlanar(const InterpolationOptions & io); @@ -43,15 +41,13 @@ namespace INTERP_KERNEL // Main function to interpolate triangular and quadratic meshes template - int interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result, const char *method); - + int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method); public: bool doRotate() const { return asLeafInterpPlanar().doRotate(); } double medianPlane() const { return asLeafInterpPlanar().medianPlane(); } template void performAdjustmentOfBB(PlanarIntersector* intersector, std::vector& bbox) const { return asLeafInterpPlanar().performAdjustmentOfBB(intersector,bbox); } - protected: RealPlanar& asLeafInterpPlanar() { return static_cast(*this); } const RealPlanar& asLeafInterpPlanar() const { return static_cast< const RealPlanar& >(*this); } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index 1e720fd25..49eed13aa 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -20,6 +20,7 @@ #define __INTERPOLATIONPLANAR_TXX__ #include "InterpolationPlanar.hxx" +#include "Interpolation.txx" #include "InterpolationOptions.hxx" #include "PlanarIntersector.hxx" #include "PlanarIntersector.txx" @@ -36,10 +37,6 @@ namespace INTERP_KERNEL { - - template - const double InterpolationPlanar::DEFAULT_PRECISION=1.e-12; - /** * \defgroup interpolationPlanar InterpolationPlanar * @@ -59,7 +56,6 @@ namespace INTERP_KERNEL { } - /** * \brief Function used to set the options for the intersection calculation * \details The following options can be modified: diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 0403e79f9..042c2c49a 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -41,6 +41,7 @@ Geometric2DIntersector.txx \ INTERPKERNELDefines.hxx \ InterpKernelMatrix.hxx \ Interpolation.hxx \ +Interpolation.txx \ Interpolation2D.hxx \ Interpolation2D.txx \ Interpolation3D.hxx \ @@ -91,6 +92,8 @@ TransformedTriangleInline.hxx \ TranslationRotationMatrix.hxx \ TriangulationIntersector.hxx \ TriangulationIntersector.txx \ +IntegralUniformIntersector.hxx \ +IntegralUniformIntersector.txx \ UnitTetraIntersectionBary.hxx \ VTKNormalizedUnstructuredMesh.hxx \ VTKNormalizedUnstructuredMesh.txx \ @@ -104,7 +107,9 @@ PlanarIntersectorP0P1.hxx \ PlanarIntersectorP0P1.txx \ PlanarIntersectorP1P0.hxx \ PlanarIntersectorP1P0.txx \ -VolSurfFormulae.hxx +VolSurfFormulae.hxx \ +VolSurfUser.hxx \ +VolSurfUser.txx # Libraries targets diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx index 52c49856a..7effb8eff 100644 --- a/src/INTERP_KERNEL/VolSurfFormulae.hxx +++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx @@ -16,8 +16,8 @@ // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // -#ifndef VOLSURFFORMULAE -#define VOLSURFFORMULAE +#ifndef __VOLSURFFORMULAE_HXX__ +#define __VOLSURFFORMULAE_HXX__ #include diff --git a/src/INTERP_KERNEL/VolSurfUser.hxx b/src/INTERP_KERNEL/VolSurfUser.hxx new file mode 100644 index 000000000..ef89fcd53 --- /dev/null +++ b/src/INTERP_KERNEL/VolSurfUser.hxx @@ -0,0 +1,33 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __VOLSURFUSER_HXX__ +#define __VOLSURFUSER_HXX__ + +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + template + double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords); + + template + double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim); +} + +#endif diff --git a/src/INTERP_KERNEL/VolSurfUser.txx b/src/INTERP_KERNEL/VolSurfUser.txx new file mode 100644 index 000000000..e65779601 --- /dev/null +++ b/src/INTERP_KERNEL/VolSurfUser.txx @@ -0,0 +1,168 @@ +// Copyright (C) 2007-2008 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. +// +// 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 +// +#ifndef __VOLSURFUSER_TXX__ +#define __VOLSURFUSER_TXX__ + +#include "VolSurfUser.hxx" +#include "VolSurfFormulae.hxx" +#include "InterpolationUtils.hxx" + +namespace INTERP_KERNEL +{ + template + double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords) + { + switch(type) + { + case INTERP_KERNEL::NORM_TRI3 : + case INTERP_KERNEL::NORM_TRI6 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + + return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1), + coords+(SPACEDIM*N2), + coords+(SPACEDIM*N3), + SPACEDIM); + } + break; + + case INTERP_KERNEL::NORM_QUAD4 : + case INTERP_KERNEL::NORM_QUAD8 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + int N4 = OTT::coo2C(connec[3]); + + return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1, + coords+SPACEDIM*N2, + coords+SPACEDIM*N3, + coords+SPACEDIM*N4, + SPACEDIM); + } + break; + + case INTERP_KERNEL::NORM_POLYGON : + { + const double **pts=new const double *[lgth]; + for(int inod=0;inod::coo2C(connec[inod]); + double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM); + delete [] pts; + return val; + } + break; + case INTERP_KERNEL::NORM_TETRA4 : + case INTERP_KERNEL::NORM_TETRA10 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + int N4 = OTT::coo2C(connec[3]); + + return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1, + coords+SPACEDIM*N2, + coords+SPACEDIM*N3, + coords+SPACEDIM*N4); + } + break; + + case INTERP_KERNEL::NORM_PYRA5 : + case INTERP_KERNEL::NORM_PYRA13 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + int N4 = OTT::coo2C(connec[3]); + int N5 = OTT::coo2C(connec[4]); + + return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1, + coords+SPACEDIM*N2, + coords+SPACEDIM*N3, + coords+SPACEDIM*N4, + coords+SPACEDIM*N5); + } + break; + + case INTERP_KERNEL::NORM_PENTA6 : + case INTERP_KERNEL::NORM_PENTA15 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + int N4 = OTT::coo2C(connec[3]); + int N5 = OTT::coo2C(connec[4]); + int N6 = OTT::coo2C(connec[5]); + + return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1, + coords+SPACEDIM*N2, + coords+SPACEDIM*N3, + coords+SPACEDIM*N4, + coords+SPACEDIM*N5, + coords+SPACEDIM*N6); + } + break; + + case INTERP_KERNEL::NORM_HEXA8 : + case INTERP_KERNEL::NORM_HEXA20 : + { + int N1 = OTT::coo2C(connec[0]); + int N2 = OTT::coo2C(connec[1]); + int N3 = OTT::coo2C(connec[2]); + int N4 = OTT::coo2C(connec[3]); + int N5 = OTT::coo2C(connec[4]); + int N6 = OTT::coo2C(connec[5]); + int N7 = OTT::coo2C(connec[6]); + int N8 = OTT::coo2C(connec[7]); + + return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1, + coords+SPACEDIM*N2, + coords+SPACEDIM*N3, + coords+SPACEDIM*N4, + coords+SPACEDIM*N5, + coords+SPACEDIM*N6, + coords+SPACEDIM*N7, + coords+SPACEDIM*N8); + } + break; + + case INTERP_KERNEL::NORM_POLYHED : + { + throw INTERP_KERNEL::Exception("Polyedra Not yet implemented !"); + } + break; + default: + throw INTERP_KERNEL::Exception("Not recognized cell type to get Area/Volume on it !"); + } + } + + template + double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim) + { + if(spaceDim==3) + return computeVolSurfOfCell(type,connec,lgth,coords); + if(spaceDim==2) + return computeVolSurfOfCell(type,connec,lgth,coords); + throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 2 or 3"); + } +} + +#endif diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx index 33f7214c0..19deae8d1 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.cxx +++ b/src/MEDCoupling/MEDCouplingCMesh.cxx @@ -175,7 +175,7 @@ void MEDCouplingCMesh::getBoundingBox(double *bbox) const //not implemented yet ! } -MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField() const +MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const { //not implemented yet ! return 0; diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx index 212535c44..59fc9c2e7 100644 --- a/src/MEDCoupling/MEDCouplingCMesh.hxx +++ b/src/MEDCoupling/MEDCouplingCMesh.hxx @@ -45,7 +45,7 @@ namespace ParaMEDMEM DataArrayDouble *coordsZ=0); // tools void getBoundingBox(double *bbox) const; - MEDCouplingFieldDouble *getMeasureField() const; + MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; private: MEDCouplingCMesh(); ~MEDCouplingCMesh(); diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 5b3e5cac6..d1c82fc26 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -86,9 +86,9 @@ void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMe } } -MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh) const +MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const { - return mesh->getMeasureField(); + return mesh->getMeasureField(isAbs); } /*! @@ -149,7 +149,7 @@ void MEDCouplingFieldDiscretizationP1::checkCoherencyBetween(const MEDCouplingMe } } -MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh) const +MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const { //not implemented yet. //Dual mesh to build diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index 6061bf20f..f035ba3ed 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -42,7 +42,7 @@ namespace ParaMEDMEM virtual int getNumberOfTuples(const MEDCouplingMesh *mesh) const = 0; virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0; virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0; - virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0; + virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0; virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0; }; @@ -56,7 +56,7 @@ namespace ParaMEDMEM int getNumberOfTuples(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); - MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; public: static const char REPR[]; @@ -73,7 +73,7 @@ namespace ParaMEDMEM int getNumberOfTuples(const MEDCouplingMesh *mesh) const; void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception); void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception); - MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const; + MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const; MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const; static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di); public: diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 6b1fa4e64..5c673c463 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -89,11 +89,11 @@ double MEDCouplingFieldDouble::accumulate(int compId) const return ret; } -double MEDCouplingFieldDouble::measureAccumulate(int compId) const +double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const { if(!_mesh) throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate"); - MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh); + MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh,isWAbs); const double *ptr=weight->getArray()->getConstPointer(); int nbOfValues=weight->getArray()->getNbOfElems(); double ret=0.; diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index a05035b0f..0f37543af 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -43,7 +43,7 @@ namespace ParaMEDMEM void setArray(DataArrayDouble *array); DataArrayDouble *getArray() const { return _time_discr->getArray(); } double accumulate(int compId) const; - double measureAccumulate(int compId) const; + double measureAccumulate(int compId, bool isWAbs) const; void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception); void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception); //! \b temporary diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx index 61324c2f8..597088efe 100644 --- a/src/MEDCoupling/MEDCouplingMesh.hxx +++ b/src/MEDCoupling/MEDCouplingMesh.hxx @@ -50,7 +50,7 @@ namespace ParaMEDMEM virtual int getMeshDimension() const = 0; // tools virtual void getBoundingBox(double *bbox) const = 0; - virtual MEDCouplingFieldDouble *getMeasureField() const = 0; + virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0; protected: MEDCouplingMesh() { } MEDCouplingMesh(const MEDCouplingMesh& other):_name(other._name) { } diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 4fa962cb2..a560bef91 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -19,7 +19,7 @@ #include "MEDCouplingUMesh.hxx" #include "MEDCouplingFieldDouble.hxx" #include "CellModel.hxx" -#include "VolSurfFormulae.hxx" +#include "VolSurfUser.txx" #include #include @@ -755,185 +755,29 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start * param field field on which cells the volumes are required * return field containing the volumes, area or length depending the meshdimension. */ -MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const -{ - int ipt, type ; - int nbelem = getNumberOfCells() ; - int dim_mesh = getMeshDimension(); - int dim_space = getSpaceDimension() ; - const double *coords = getCoords()->getConstPointer() ; - const int *connec = getNodalConnectivity()->getConstPointer() ; - const int *connec_index = getNodalConnectivityIndex()->getConstPointer() ; - - - MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS); - DataArrayDouble* array = DataArrayDouble::New() ; - array->alloc(nbelem, 1) ; - double *area_vol = array->getPointer() ; - - switch (dim_mesh) +MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const +{ + int ipt; + INTERP_KERNEL::NormalizedCellType type; + int nbelem=getNumberOfCells(); + int dim_space=getSpaceDimension(); + const double *coords=getCoords()->getConstPointer(); + const int *connec=getNodalConnectivity()->getConstPointer(); + const int *connec_index=getNodalConnectivityIndex()->getConstPointer(); + + MEDCouplingFieldDouble* field=MEDCouplingFieldDouble::New(ON_CELLS); + DataArrayDouble* array=DataArrayDouble::New(); + array->alloc(nbelem,1); + double *area_vol = array->getPointer(); + for(int iel=0;iel(type,connec+ipt+1,connec_index[iel+1]-ipt-1,coords,dim_space); } - + if(isAbs) + for(int iel=0;ielsetArray(array) ; array->decrRef(); field->setMesh(const_cast(this)); diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 160a57ddb..66000673d 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -67,7 +67,7 @@ namespace ParaMEDMEM MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; void renumberConnectivity(const int *newNodeNumbers); void giveElemsInBoundingBox(const double *bbox, double eps, std::vector& elems); - MEDCouplingFieldDouble *getMeasureField() const; + MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; private: MEDCouplingUMesh(); MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy); diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index 19fa0abd5..fc20b2747 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -242,7 +242,7 @@ void MEDCouplingUMeshDesc::renumberConnectivity(const int *newNodeNumbers) { } -MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const +MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField(bool isAbs) const { //not implemented yet. return 0; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index b7e672435..ca3c9b6d4 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -53,7 +53,7 @@ namespace ParaMEDMEM void findBoundaryNodes(std::vector& nodes) const; MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const; void renumberConnectivity(const int *newNodeNumbers); - MEDCouplingFieldDouble *getMeasureField() const; + MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; DataArrayInt *zipCoordsTraducer(); private: MEDCouplingUMeshDesc(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx index c00489298..ec673e640 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx @@ -1479,6 +1479,177 @@ void MEDCouplingBasicsTest::test3DInterpP0P0Empty() targetMesh->decrRef(); } +void MEDCouplingBasicsTest::test2DInterpP0IntegralUniform() +{ + MEDCouplingUMesh *targetMesh=build2DTargetMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation2D myInterpolator; + vector > res; + CPPUNIT_ASSERT_EQUAL(5,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + res.clear(); + CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[3][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[4][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + res.clear(); + targetMesh->decrRef(); + // + targetMesh=build2DTargetMeshPerm_1(); + MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper2(targetMesh); + INTERP_KERNEL::Interpolation2D myInterpolator2; + CPPUNIT_ASSERT(myInterpolator2.getMeasureAbsStatus()); + CPPUNIT_ASSERT_EQUAL(5,myInterpolator2.fromIntegralUniform(targetWrapper2,res,"P0")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12); + res.clear(); + myInterpolator2.setMeasureAbsStatus(false); + CPPUNIT_ASSERT(!myInterpolator2.getMeasureAbsStatus()); + CPPUNIT_ASSERT_EQUAL(5,myInterpolator2.fromIntegralUniform(targetWrapper2,res,"P0")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125,res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.75,sumAll(res),1e-12); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test3DSurfInterpP0IntegralUniform() +{ + MEDCouplingUMesh *targetMesh=build3DSurfTargetMesh_1(); + INTERP_KERNEL::Interpolation3DSurf myInterpolator; + MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh); + vector > res; + CPPUNIT_ASSERT_EQUAL(5,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][3],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][4],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12); + res.clear(); + CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(5,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[4][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test3DInterpP0IntegralUniform() +{ + MEDCouplingUMesh *targetMesh=build3DTargetMesh_1(); + INTERP_KERNEL::Interpolation3D myInterpolator; + MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh); + vector > res; + CPPUNIT_ASSERT_EQUAL(8,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000.,res[0][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][1],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][2],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][3],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][4],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][5],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][6],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3375000.,res[0][7],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6); + res.clear(); + CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0")); + CPPUNIT_ASSERT_EQUAL(8,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(125000.,res[0][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[1][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[2][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[3][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[4][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[5][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[6][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(3375000.,res[7][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6); + res.clear(); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test2DInterpP1IntegralUniform() +{ + MEDCouplingUMesh *targetMesh=build2DSourceMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh); + INTERP_KERNEL::Interpolation2D myInterpolator; + vector > res; + CPPUNIT_ASSERT_EQUAL(4,myInterpolator.fromIntegralUniform(targetWrapper,res,"P1")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[0][1],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[0][2],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][3],1e-12); + res.clear(); + CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P1")); + CPPUNIT_ASSERT_EQUAL(4,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[1][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[2][0],1e-12); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[3][0],1e-12); + res.clear(); + targetMesh->decrRef(); +} + +void MEDCouplingBasicsTest::test3DInterpP1IntegralUniform() +{ + MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1(); + // + MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(sourceMesh); + INTERP_KERNEL::Interpolation3D myInterpolator; + vector > res; + CPPUNIT_ASSERT_EQUAL(9,myInterpolator.fromIntegralUniform(targetWrapper,res,"P1")); + CPPUNIT_ASSERT_EQUAL(1,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][1],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[0][2],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][3],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][4],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[0][5],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][6],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][7],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2000000.,res[0][8],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6); + res.clear(); + CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P1")); + CPPUNIT_ASSERT_EQUAL(9,(int)res.size()); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[1][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[2][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[3][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[4][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[5][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[6][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[7][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(2000000.,res[8][0],1e-6); + CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6); + sourceMesh->decrRef(); +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build2DSourceMesh_1() { double sourceCoords[8]={-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7}; @@ -1517,6 +1688,27 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_1() return targetMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMeshPerm_1() +{ + double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 }; + int targetConn[18]={0,3,4,1, 1,2,4, 4,5,2, 6,7,4,3, 7,8,5,4}; + MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New(); + targetMesh->setMeshDimension(2); + targetMesh->allocateCells(5); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+4); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+7); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+10); + targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+14); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,2); + std::copy(targetCoords,targetCoords+18,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + return targetMesh; +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_2() { double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 }; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 49ddbc162..03941d922 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -63,6 +63,11 @@ namespace ParaMEDMEM CPPUNIT_TEST( test3DInterpP0P1_1 ); CPPUNIT_TEST( test3DInterpP1P0_1 ); CPPUNIT_TEST( test3DInterpP0P0Empty ); + CPPUNIT_TEST( test2DInterpP0IntegralUniform ); + CPPUNIT_TEST( test3DSurfInterpP0IntegralUniform ); + CPPUNIT_TEST( test3DInterpP0IntegralUniform ); + CPPUNIT_TEST( test2DInterpP1IntegralUniform ); + CPPUNIT_TEST( test3DInterpP1IntegralUniform ); CPPUNIT_TEST_SUITE_END(); public: void testArray(); @@ -97,9 +102,15 @@ namespace ParaMEDMEM void test3DInterpP0P1_1(); void test3DInterpP1P0_1(); void test3DInterpP0P0Empty(); + void test2DInterpP0IntegralUniform(); + void test3DSurfInterpP0IntegralUniform(); + void test3DInterpP0IntegralUniform(); + void test2DInterpP1IntegralUniform(); + void test3DInterpP1IntegralUniform(); private: MEDCouplingUMesh *build2DSourceMesh_1(); MEDCouplingUMesh *build2DTargetMesh_1(); + MEDCouplingUMesh *build2DTargetMeshPerm_1(); MEDCouplingUMesh *build2DTargetMesh_2(); MEDCouplingUMesh *build3DSurfSourceMesh_1(); MEDCouplingUMesh *build3DSurfSourceMesh_2(); diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx index c564438bb..8bb1f2d27 100644 --- a/src/ParaMEDMEM/DEC.cxx +++ b/src/ParaMEDMEM/DEC.cxx @@ -188,12 +188,12 @@ namespace ParaMEDMEM \f] */ - void DEC::renormalizeTargetField() + void DEC::renormalizeTargetField(bool isWAbs) { if (_source_group->containsMyRank()) for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++) { - double total_norm = _local_field->getVolumeIntegral(icomp+1); + double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs); double source_norm = total_norm; _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast(_union_group)->getComm()); @@ -202,7 +202,7 @@ namespace ParaMEDMEM { for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++) { - double total_norm = _local_field->getVolumeIntegral(icomp+1); + double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs); double source_norm=total_norm; _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast(_union_group)->getComm()); diff --git a/src/ParaMEDMEM/DEC.hxx b/src/ParaMEDMEM/DEC.hxx index 0f4c545bb..d6cc17a7b 100644 --- a/src/ParaMEDMEM/DEC.hxx +++ b/src/ParaMEDMEM/DEC.hxx @@ -50,7 +50,7 @@ namespace ParaMEDMEM virtual void synchronize()=0; virtual ~DEC(); virtual void computeProcGroup() { } - void renormalizeTargetField(); + void renormalizeTargetField(bool isWAbs); protected: void compareFieldAndMethod() const throw(INTERP_KERNEL::Exception); protected: diff --git a/src/ParaMEDMEM/InterpKernelDEC.cxx b/src/ParaMEDMEM/InterpKernelDEC.cxx index 6badf5117..444378a1f 100644 --- a/src/ParaMEDMEM/InterpKernelDEC.cxx +++ b/src/ParaMEDMEM/InterpKernelDEC.cxx @@ -222,7 +222,7 @@ namespace ParaMEDMEM { _interpolation_matrix->multiply(*_local_field->getField()); if (getForcedRenormalization()) - renormalizeTargetField(); + renormalizeTargetField(getMeasureAbsStatus()); } } @@ -249,7 +249,7 @@ namespace ParaMEDMEM _interpolation_matrix->multiply(*_local_field->getField()); if (getForcedRenormalization()) - renormalizeTargetField(); + renormalizeTargetField(getMeasureAbsStatus()); } else if (_target_group->containsMyRank()) diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx index 6edb3d40e..42b6e8d63 100644 --- a/src/ParaMEDMEM/InterpolationMatrix.cxx +++ b/src/ParaMEDMEM/InterpolationMatrix.cxx @@ -154,7 +154,7 @@ namespace ParaMEDMEM MEDCouplingFieldDouble *target_triangle_surf=0; if(needTargetSurf) - target_triangle_surf = distant_support.getMeasureField(); + target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus()); fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf); if(needTargetSurf) @@ -320,7 +320,7 @@ namespace ParaMEDMEM void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator) { - MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(); + MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus()); _deno_multiply.resize(_coeffs.size()); vector >::iterator iter6=_deno_multiply.begin(); const double *values=source_triangle_surf->getArray()->getConstPointer(); diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx index 6afb23911..58532760a 100644 --- a/src/ParaMEDMEM/ParaFIELD.cxx +++ b/src/ParaMEDMEM/ParaFIELD.cxx @@ -211,10 +211,10 @@ namespace ParaMEDMEM /*! This method retrieves the integral of component \a icomp over the all domain. */ - double ParaFIELD::getVolumeIntegral(int icomp) const + double ParaFIELD::getVolumeIntegral(int icomp, bool isWAbs) const { CommInterface comm_interface = _topology->getProcGroup()->getCommInterface(); - double integral=_field->measureAccumulate(icomp); + double integral=_field->measureAccumulate(icomp,isWAbs); double total=0.; const MPI_Comm* comm = (dynamic_cast(_topology->getProcGroup()))->getComm(); comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm); diff --git a/src/ParaMEDMEM/ParaFIELD.hxx b/src/ParaMEDMEM/ParaFIELD.hxx index a9bcf57b8..22adbf5c0 100644 --- a/src/ParaMEDMEM/ParaFIELD.hxx +++ b/src/ParaMEDMEM/ParaFIELD.hxx @@ -48,7 +48,7 @@ namespace ParaMEDMEM Topology* getTopology() const { return _topology; } ParaMESH* getSupport() const { return _support; } int nbComponents() const; - double getVolumeIntegral(int icomp) const; + double getVolumeIntegral(int icomp, bool isWAbs) const; double getL2Norm()const { return -1; } private: MEDCouplingFieldDouble* _field; diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx index cfaef00ff..ed54396b6 100644 --- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx +++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx @@ -199,7 +199,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_2D_(const char *srcMeth, const char *ta if (source_group->containsMyRank()) { - field_before_int = parafield->getVolumeIntegral(0); + field_before_int = parafield->getVolumeIntegral(0,true); dec.synchronize(); cout<<"DEC usage"<myRank()+1; aRemover.Register(filename.str().c_str()); - field_after_int = parafield->getVolumeIntegral(0); + field_after_int = parafield->getVolumeIntegral(0,true); // MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD); @@ -406,7 +406,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_3D_(const char *srcMeth, const char *ta if (source_group->containsMyRank()) { - field_before_int = parafield->getVolumeIntegral(0); + field_before_int = parafield->getVolumeIntegral(0,true); dec.synchronize(); cout<<"DEC usage"<myRank()+1; aRemover.Register(filename.str().c_str()); - field_after_int = parafield->getVolumeIntegral(0); + field_after_int = parafield->getVolumeIntegral(0,true); CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6); @@ -1190,7 +1190,7 @@ void ParaMEDMEMTest::testAsynchronousInterpKernelDEC_2D(double dtA, double tmaxA cout << "testAsynchronousInterpKernelDEC_2D" << rank << " time " << time << " dtB " << dtB << " tmaxB " << tmaxB << endl ; dec.recvData( time ); - double vi = parafield->getVolumeIntegral(0); + double vi = parafield->getVolumeIntegral(0,true); cout << "testAsynchronousInterpKernelDEC_2D" << rank << " time " << time << " VolumeIntegral " << vi << " time*10000 " << time*10000 << endl ; diff --git a/src/ParaMEDMEM/Test/test_perf.cxx b/src/ParaMEDMEM/Test/test_perf.cxx index 72ac18b65..f5872c15b 100644 --- a/src/ParaMEDMEM/Test/test_perf.cxx +++ b/src/ParaMEDMEM/Test/test_perf.cxx @@ -210,7 +210,7 @@ void testInterpKernelDEC_2D(const string& filename_xml1, const string& meshname1 double field_after_int; if (source_group->containsMyRank()){ - field_before_int = parafield->getVolumeIntegral(0); + field_before_int = parafield->getVolumeIntegral(0,true); get_time( &telps, &tcpu_u, &tcpu_s, &tcpu ); dec.synchronize(); get_time( &telps, &tcpu_u, &tcpu_s, &tcpu ); @@ -229,7 +229,7 @@ void testInterpKernelDEC_2D(const string& filename_xml1, const string& meshname1 cout << "SEND DATA : Telapse = " << telps << " TuserCPU = " << tcpu_u << " TsysCPU = " << tcpu_s << " TCPU = " << tcpu << endl; dec.recvData(); - field_after_int = parafield->getVolumeIntegral(0); + field_after_int = parafield->getVolumeIntegral(0,true); // CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, epsilon); }