From: Anthony Geay Date: Wed, 3 Jan 2018 09:57:39 +0000 (+0100) Subject: First draft OK X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=d4f3f3d409040369b1368874be891a4d5f93dec5;p=tools%2Fmedcoupling.git First draft OK --- diff --git a/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx b/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx new file mode 100644 index 000000000..51daf4284 --- /dev/null +++ b/src/INTERP_KERNEL/Bases/InterpKernelAssert.hxx @@ -0,0 +1,34 @@ +// Copyright (C) 2018 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 __INTERPKERNELASSERT_HXX__ +#define __INTERPKERNELASSERT_HXX__ + +#include "InterpKernelException.hxx" + +#include + +#define IKAssert(a) { bool verdict(a); \ + if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " !"; throw INTERP_KERNEL::Exception(osszz.str()); } } + +#define IKAssertMsg(a,b) { bool verdict(a); \ + if(!verdict) { std::ostringstream osszz; osszz << "Assertion \"" << #a << "\" failed into " << __FILE__ << " at line " << __LINE__ << " with message \"" << b << "\" !"; throw INTERP_KERNEL::Exception(osszz.str()); } } + +#endif diff --git a/src/INTERP_KERNEL/BoundingBox.cxx b/src/INTERP_KERNEL/BoundingBox.cxx index eeccb1ee0..b20ca4100 100644 --- a/src/INTERP_KERNEL/BoundingBox.cxx +++ b/src/INTERP_KERNEL/BoundingBox.cxx @@ -164,7 +164,12 @@ namespace INTERP_KERNEL void BoundingBox::toCompactData(double data[6]) const { - std::copy(_coords,_coords+6,data); + data[0]=_coords[XMIN]; + data[1]=_coords[XMAX]; + data[2]=_coords[YMIN]; + data[3]=_coords[YMAX]; + data[4]=_coords[ZMIN]; + data[5]=_coords[ZMAX]; } } diff --git a/src/INTERP_KERNEL/Interpolation1D0D.txx b/src/INTERP_KERNEL/Interpolation1D0D.txx index db6296b5b..3227d3a39 100644 --- a/src/INTERP_KERNEL/Interpolation1D0D.txx +++ b/src/INTERP_KERNEL/Interpolation1D0D.txx @@ -32,6 +32,8 @@ #include "BBTree.txx" +#include "InterpKernelAssert.hxx" + namespace INTERP_KERNEL { /** @@ -43,16 +45,17 @@ namespace INTERP_KERNEL { //#if __cplusplus >= 201103L constexpr int SPACEDIM=MyMeshType::MY_SPACEDIM; + using ConnType=typename MyMeshType::MyConnType; + IKAssert(SPACEDIM==3); if(InterpolationOptions::getIntersectionType() != PointLocator) INTERP_KERNEL::Exception("Invalid 1D/0D intersection type specified : must be PointLocator."); std::string methC ( InterpolationOptions::filterInterpolationMethod(method) ); - if(methC=="P1P1") + if(methC!="P1P1") throw Exception("Invalid method chosen must be in \"P1P1\"."); - - typedef typename MyMeshType::MyConnType ConnType; + const double epsilon(getPrecision()); // create MeshElement objects corresponding to each element of the two meshes const unsigned long numSrcElems(srcMesh.getNumberOfElements()), numTargetElems(targetMesh.getNumberOfElements()); @@ -80,32 +83,37 @@ namespace INTERP_KERNEL } adjustBoundingBoxes(bboxes); - const double *bboxPtr=0; + const double *bboxPtr(nullptr); if(numSrcElems>0) bboxPtr=&bboxes[0]; - BBTree<3,ConnType> tree(bboxPtr, srcElemIdx, 0, numSrcElems); - - // for each target element, get source elements with which to calculate intersection - // - calculate intersection by calling intersectCells + BBTree tree(bboxPtr, srcElemIdx, 0, numSrcElems); + const ConnType *trgConnPtr(targetMesh.getConnectivityPtr()),*trgConnIPtr(targetMesh.getConnectivityIndexPtr()); + const ConnType *srcConnPtr(srcMesh.getConnectivityPtr()),*srcConnIPtr(srcMesh.getConnectivityIndexPtr()); + const double *trgCooPtr(targetMesh.getCoordinatesPtr()),*srcCooPtr(srcMesh.getCoordinatesPtr()); for(unsigned long i = 0; i < numTargetElems; ++i) { - std::vector intersectElems; - - //tree.getElementsAroundPoint(targetBox, intersectElems); - - //if ( !intersectElems.empty() ) - // intersector->intersectCells(targetIdx,intersectElems,result); + IKAssert(trgConnIPtr[i+1]==i+1 && trgConnIPtr[i]==i); + std::vector srcSegCondidates; + const double *trgCellPosition(trgCooPtr+SPACEDIM*trgConnPtr[i]); + typename MatrixType::value_type& resRow(result[i]); + tree.getElementsAroundPoint(trgCellPosition, srcSegCondidates); + for(auto srcSeg: srcSegCondidates) + { + IKAssertMsg(srcConnIPtr[srcSeg+1]==2*(srcSeg+1) && srcConnIPtr[srcSeg]==2*srcSeg,"Only implemented for linear 1D source"); + double bc0(0.),bc1(0.); + ConnType srcNode0(srcConnPtr[2*srcSeg]),srcNode1(srcConnPtr[2*srcSeg+1]); + if(IsPointOn3DSeg(srcCooPtr+SPACEDIM*srcNode0,srcCooPtr+SPACEDIM*srcNode1,trgCellPosition,epsilon,bc0,bc1)) + { + resRow.insert(std::make_pair(srcNode0,bc0)); + resRow.insert(std::make_pair(srcNode1,bc1)); + continue; + } + } } - - // free allocated memory delete [] srcElemIdx; - for(unsigned long i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } + delete srcElems[i]; return srcMesh.getNumberOfNodes(); - } //#endif } diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index e222e1b07..fad2faea7 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -68,7 +68,7 @@ namespace INTERP_KERNEL /* calcul la surface d'un triangle */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ - inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3) + inline double Surf_Tri(const double *P_1,const double *P_2,const double *P_3) { double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]); double Surface = 0.5*fabs(A); @@ -83,9 +83,9 @@ namespace INTERP_KERNEL //fonction qui calcul le determinant des vecteurs: P3P1 et P3P2 //(cf doc CGAL). - inline double mon_determinant(const double* P_1, - const double* P_2, - const double* P_3) + inline double mon_determinant(const double *P_1, + const double *P_2, + const double *P_3) { double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]); return mon_det; @@ -98,8 +98,7 @@ namespace INTERP_KERNEL { double X=P_1[0]-P_2[0]; double Y=P_1[1]-P_2[1]; - double norme=sqrt(X*X+Y*Y); - return norme; + return sqrt(X*X+Y*Y); } /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -536,6 +535,29 @@ namespace INTERP_KERNEL } } + /*! + * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space. + * If point \a point is further from S than eps false is returned. + * If point \a point projection on S is outside S false is also returned. + * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set. + */ + inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1) + { + double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]}; + double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2])); + double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB)); + double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]}; + double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]}; + double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2])); + if(dist_P_AB>=eps) + return false;//to far from segment [segStart,segStop] + if(AP_dot_AB<-eps || AP_dot_AB>1.+eps) + return false; + AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.); + bc0=1.-AP_dot_AB; bc1=AP_dot_AB; + return true; + } + /*! * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices. * This method makes the assumption that: diff --git a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx index ce207679e..111b6186c 100644 --- a/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx +++ b/src/MEDCoupling/MEDCouplingNormalizedUnstructuredMesh.txx @@ -22,6 +22,7 @@ #define __MEDCOUPLINGNORMALIZEDUNSTRUCTUREDMESH_TXX__ #include "MEDCouplingNormalizedUnstructuredMesh.hxx" +#include "InterpKernelAssert.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCoupling1GTUMesh.hxx" @@ -124,6 +125,7 @@ MEDCouplingNormalizedUnstructuredMesh::~MEDCouplingNormalizedU template void MEDCouplingNormalizedUnstructuredMesh::prepare() { + IKAssert(_mesh->getSpaceDimension()==SPACEDIM); const MEDCoupling::MEDCouplingUMesh *m1(dynamic_cast(_mesh)); if(m1) { diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx index 71d30669a..683a8d195 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -1465,7 +1465,7 @@ void MEDCouplingBasicsTestInterp::test3DInterpP0P0Empty() sourceMesh->allocateCells(0); sourceMesh->finishInsertingCells(); DataArrayDouble *myCoords=DataArrayDouble::New(); - myCoords->alloc(0,0); + myCoords->alloc(0,2); sourceMesh->setCoords(myCoords); myCoords->decrRef(); MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();