From: abn Date: Mon, 18 Feb 2019 10:29:55 +0000 (+0100) Subject: Remapper: bug fix PointLocator for srcMeshDim=2, trgtMeshDim=3, spaceDim=3 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=90d41e009fcb2989fe47afa5a0493b475d6197c9;p=tools%2Fmedcoupling.git Remapper: bug fix PointLocator for srcMeshDim=2, trgtMeshDim=3, spaceDim=3 This is now consistent with the 1D/2D situation: we revert the barycenter localisation. Source mesh barycenters are localised in target mesh. --- diff --git a/src/INTERP_KERNEL/Interpolation2D3D.txx b/src/INTERP_KERNEL/Interpolation2D3D.txx index c9c7dad82..a9b725243 100644 --- a/src/INTERP_KERNEL/Interpolation2D3D.txx +++ b/src/INTERP_KERNEL/Interpolation2D3D.txx @@ -99,11 +99,8 @@ namespace INTERP_KERNEL intersectFaces, getSplittingPolicy()); break; - case PointLocator:// switch target and source - intersector=new PointLocator3DIntersectorP0P0(srcMesh,targetMesh,getPrecision()); - break; default: - throw INTERP_KERNEL::Exception("Invalid 3D to 2D intersection type for P0P0 interp specified : must be Triangulation or PointLocator."); + throw INTERP_KERNEL::Exception("Invalid 2D to 3D intersection type for P0P0 interp specified : must be Triangulation."); } } else diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx index bec0b62e7..f914b5e40 100644 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -552,19 +552,33 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() } else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3) { - MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh); - MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); - INTERP_KERNEL::Interpolation2D3D interpolation(*this); - nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); - INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); - if(!duplicateFaces.empty()) + if(getIntersectionType()==INTERP_KERNEL::PointLocator) + { + MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); + INTERP_KERNEL::Interpolation3D interpolation(*this); + std::vector > matrixTmp; + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); + ReverseMatrix(matrixTmp,nbCols,_matrix); + nbCols=matrixTmp.size(); + } + else { - std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n"; - for(std::map >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++) + MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh); + MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); + INTERP_KERNEL::Interpolation2D3D interpolation(*this); + nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method); + INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); + if(!duplicateFaces.empty()) { - oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : "; - std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator(oss," ")); - oss << std::endl; + std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n"; + for(std::map >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++) + { + oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : "; + std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator(oss," ")); + oss << std::endl; + } } } } diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py index a2efa36ab..3d05957ca 100644 --- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py @@ -988,18 +988,34 @@ class MEDCouplingBasicsTest(unittest.TestCase): # mt=MEDCouplingUMesh("target",2) ; mt.allocateCells() mt.insertNextCell(NORM_TRI3,[0,4,6]) - mt.insertNextCell(NORM_TRI3,[1,5,7]) mt.setCoords(ms.getCoords()[:]) mt.zipCoords() - # + # 3D to 2D rem=MEDCouplingRemapper() rem.setIntersectionType(PointLocator) rem.prepare(ms,mt,"P0P0") - self.assertEqual(rem.getCrudeMatrix(),[{0: 1.0}, {1: 1.0}]) + self.assertEqual(rem.getCrudeMatrix(),[{0: 1.0}]) + + mt1=MEDCouplingUMesh("target",2) ; mt1.allocateCells() + mt1.insertNextCell(NORM_TRI3,[0,4,6]) + mt1.insertNextCell(NORM_TRI3,[1,5,7]) + mt1.setCoords(ms.getCoords()[:]) + mt1.zipCoords() + rem1=MEDCouplingRemapper() + rem1.setIntersectionType(PointLocator) + rem1.prepare(ms,mt1,"P0P0") + self.assertEqual(rem1.getCrudeMatrix(),[{0: 1.0}, {1: 1.0}]) + + # 2D to 3D rem2=MEDCouplingRemapper() rem2.setIntersectionType(PointLocator) - rem2.prepare(mt,ms,"P0P0") # reverse mt<->ms - self.assertEqual(rem2.getCrudeMatrix(),[{0: 1.0}, {1: 1.0}]) + rem2.prepare(mt,ms,"P0P0") # was segfaulting + self.assertEqual(rem2.getCrudeMatrix(),[{0: 1.0}, {}]) + + rem3=MEDCouplingRemapper() + rem3.setIntersectionType(PointLocator) + rem3.prepare(mt1,ms,"P0P0") + self.assertEqual(rem3.getCrudeMatrix(),[{0: 1.0}, {1: 1.0}]) pass def test2D1Dand1D2DPointLocator1(self): @@ -1192,7 +1208,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): rem.setIntersectionType(PointLocator) self.assertEqual(rem.prepare(src,trg,"P1P1"),1) mat=rem.getCrudeCSRMatrix() - row=array([2,2, 3,3, 0,0, 1,1]) # here ref to target point 3 + row=array([2,2, 3,3, 0,0, 1,1]) # here ref to target point 3 col=array([1,2, 1,2, 1,2, 1,2]) data=array([0.1,0.9, 0.3,0.7, 0.5,0.5, 0.8,0.2]) mExp2=csr_matrix((data,(row,col)),shape=(4,3)) @@ -1299,7 +1315,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): ref=float(m.getMeasureField(True).getArray()) self.assertTrue(abs(res-ref)/ref<1e-12) pass - + def checkMatrix(self,mat1,mat2,nbCols,eps): self.assertEqual(len(mat1),len(mat2)) for i in range(len(mat1)):