From: Adrien Bruneton Date: Wed, 19 Feb 2014 14:15:47 +0000 (+0100) Subject: Bug fix: bounding box for quadratic elements spaceDim=2/meshDim=1 (i.e. SEG3) X-Git-Tag: V7_4_0a1~40^2~1 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=be4c3bb042d5426fbbe54378b9d7b35173ab27ef;p=tools%2Fmedcoupling.git Bug fix: bounding box for quadratic elements spaceDim=2/meshDim=1 (i.e. SEG3) was not properly computed. Was leading to spurious bugs in Intersect2DMeshes(). --- diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index e4ecd62bc..431ee8ad9 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -70,36 +70,41 @@ QuadraticPolygon::~QuadraticPolygon() { } -QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector& nodes) +QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector& nodes, bool isClosed) { QuadraticPolygon *ret=new QuadraticPolygon; std::size_t size=nodes.size(); - for(std::size_t i=0;ipushBack(new EdgeLin(nodes[i],nodes[(i+1)%size])); nodes[i]->decrRef(); } + if(!isClosed) + nodes[size-1]->decrRef(); return ret; } -QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector& nodes) +QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector& nodes, bool isClosed) { QuadraticPolygon *ret=new QuadraticPolygon; std::size_t size=nodes.size(); - for(std::size_t i=0;ipushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)])); + ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%quad_offset])); else - ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)])); - nodes[i]->decrRef(); nodes[i+size/2]->decrRef(); + ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+quad_offset],nodes[(i+1)%quad_offset])); + nodes[i]->decrRef(); nodes[i+quad_offset]->decrRef(); } + if(!isClosed) + nodes[quad_offset-1]->decrRef(); return ret; } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index 4df9e861b..a8097a256 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -46,8 +46,8 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT QuadraticPolygon() { } INTERPKERNEL_EXPORT QuadraticPolygon(const QuadraticPolygon& other):ComposedEdge(other) { } INTERPKERNEL_EXPORT QuadraticPolygon(const char *fileName); - INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector& nodes); - INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector& nodes); + INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector& nodes, bool isClosed=true); + INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector& nodes, bool isClosed=true); INTERPKERNEL_EXPORT static void BuildDbgFile(const std::vector& nodes, const char *fileName); INTERPKERNEL_EXPORT ~QuadraticPolygon(); INTERPKERNEL_EXPORT void closeMe() const; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 29d78628e..7f339ddea 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -6272,7 +6272,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) const { int mDim(getMeshDimension()),sDim(getSpaceDimension()); - if(mDim!=2 || sDim!=2) + if(sDim!=2) // Compute refined boundary box for quadratic elements only in 2D. return getBoundingBoxForBBTreeFast(); else { @@ -6339,8 +6339,10 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const } /*! - * This method aggregate the bbox regarding foreach 2D cell in \a this the whole shape. So this method is particulary useful for 2D meshes having quadratic cells - * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes. + * This method aggregates the bbox of each 2D cell in \a this considering the whole shape. This method is particularly + * useful for 2D meshes having quadratic cells + * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers + * the two extremities of the arc of circle). * * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12) * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components. @@ -6352,7 +6354,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc { checkFullyDefined(); int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes()); - if(mDim!=2 || spaceDim!=2) + if(spaceDim!=2) throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!"); MEDCouplingAutoRefCountObjectPtr ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim); double *bbox(ret->getPointer()); @@ -6371,9 +6373,9 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]); } if(!cm.isQuadratic()) - pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes); + pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes, mDim==2); else - pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes); + pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes, mDim==2); INTERP_KERNEL::Bounds b; pol->fillBounds(b); delete pol; bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); } diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx index fb039d015..34846350c 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx @@ -1798,6 +1798,64 @@ void MEDCouplingBasicsTest5::testIntersect2DMeshesTmp6() d2->decrRef(); } +void MEDCouplingBasicsTest5::testIntersect2DMeshesTmp7() +{ + double eps = 1.0e-8; + // coordinates circle - SEE getCircle() on the Python side + DataArrayDouble *coords1=DataArrayDouble::New(); + const double coordsData1[16]={0.5328427124746189, -0.08284271247461905, -0.03284271247461901, 0.4828427124746191, -0.03284271247461906, -0.082842712474619, 0.5328427124746191, 0.482842712474619}; + coords1->useArray(coordsData1,false,CPP_DEALLOC,8,2); + // connectivity + DataArrayInt *conn1=DataArrayInt::New(); + const int connData1[5]={INTERP_KERNEL::NORM_QPOLYG,0,1,2,3}; + conn1->useArray(connData1,false,CPP_DEALLOC,5,1); + DataArrayInt *connI1=DataArrayInt::New(); + const int connIData1[2]={0,5}; + connI1->useArray(connIData1,false,CPP_DEALLOC,2,1); + MEDCouplingUMesh *m1=MEDCouplingUMesh::New("circle",2); + m1->setCoords(coords1); + m1->setConnectivity(conn1,connI1,true); + coords1->decrRef(); conn1->decrRef(); connI1->decrRef(); + + // square + DataArrayDouble *coords2=DataArrayDouble::New(); + const double coordsData2[8]={-0.5,-0.5, -0.5, 0.5, 0.5, 0.5, 0.5,-0.5}; + coords2->useArray(coordsData2,false,CPP_DEALLOC,4,2); + // connectivity + DataArrayInt *conn2=DataArrayInt::New(); + const int connData2[5]={INTERP_KERNEL::NORM_POLYGON, 0,1,2,3}; + conn2->useArray(connData2,false,CPP_DEALLOC,5,1); + DataArrayInt *connI2=DataArrayInt::New(); + const int connIData2[2]={0,5}; + connI2->useArray(connIData2,false,CPP_DEALLOC,2,1); + MEDCouplingUMesh *m2=MEDCouplingUMesh::New("square",2); + m2->setCoords(coords2); + m2->setConnectivity(conn2,connI2,true); + coords2->decrRef(); conn2->decrRef(); connI2->decrRef(); + + DataArrayInt * resToM1 = 0, * resToM2 = 0; + MEDCouplingUMesh *m_intersec=MEDCouplingUMesh::Intersect2DMeshes(m2, m1, eps, resToM1, resToM2); + m_intersec->zipCoords(); + + const double coo_tgt[34]={-0.5, -0.5, -0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.03284271247461901, 0.4828427124746191, \ + -0.014575131106459124, 0.5000000000000001, 0.5, -0.11224989991991996, 0.24271243444677046, 0.5, 0.5, 0.19387505004004, \ + -0.04799910280454185, -0.06682678787499614, -0.023843325638122054, 0.4915644577163915, 0.5, -0.30612494995996, 0.0, -0.5,\ + -0.5, 0.0, -0.25728756555322957, 0.5, -0.023843325638122026, 0.49156445771639157, -0.04799910280454181, -0.06682678787499613}; + const int conn_tgt[22]={32, 5, 2, 6, 4, 7, 8, 9, 10, 32, 6, 3, 0, 1, 5, 4, 11, 12, 13, 14, 15, 16}; + const int connI_tgt[3]={0, 9, 22}; + const int res1_tgt[2] = {0, 0}; + const int res2_tgt[2] = {0, -1}; + + CPPUNIT_ASSERT(std::equal(conn_tgt,conn_tgt+22,m_intersec->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(connI_tgt,connI_tgt+3,m_intersec->getNodalConnectivityIndex()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(res1_tgt,res1_tgt+2,resToM1->getConstPointer())); + CPPUNIT_ASSERT(std::equal(res2_tgt,res2_tgt+2,resToM2->getConstPointer())); + for(int i=0;i<34;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(coo_tgt[i],m_intersec->getCoords()->getIJ(0,i),1e-12); + m1->decrRef(); m2->decrRef(); m_intersec->decrRef(); + resToM1->decrRef(); resToM2->decrRef(); +} + void MEDCouplingBasicsTest5::testDAIBuildSubstractionOptimized1() { const int tab1[7]={1,3,5,6,7,9,13}; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx index 9df6e65fb..4b1770100 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx @@ -72,6 +72,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( testDAIPartitionByDifferentValues1 ); CPPUNIT_TEST( testDAICheckMonotonic1 ); CPPUNIT_TEST( testIntersect2DMeshesTmp6 ); + CPPUNIT_TEST( testIntersect2DMeshesTmp7 ); CPPUNIT_TEST( testDAIBuildSubstractionOptimized1 ); CPPUNIT_TEST( testDAIIsStrictlyMonotonic1 ); CPPUNIT_TEST( testSimplexize3 ); @@ -113,6 +114,7 @@ namespace ParaMEDMEM void testDAIPartitionByDifferentValues1(); void testDAICheckMonotonic1(); void testIntersect2DMeshesTmp6(); + void testIntersect2DMeshesTmp7(); void testDAIBuildSubstractionOptimized1(); void testDAIIsStrictlyMonotonic1(); void testSimplexize3(); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index 734391fa9..22b18a2f3 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -10419,6 +10419,35 @@ class MEDCouplingBasicsTest(unittest.TestCase): # Check coordinates: self.assertTrue(m3.getCoords().isEqual(m5.getCoords(), eps)) + def testIntersect2DMeshesTmp7(self): + eps = 1.0e-8 + coords = [-0.5,-0.5, -0.5, 0.5, 0.5, 0.5, 0.5,-0.5] + connec = range(4) + m1 = MEDCouplingUMesh.New("box", 2) + m1.allocateCells(1) + meshCoords = DataArrayDouble.New(coords, len(coords)/2, 2) + m1.setCoords(meshCoords) + m1.insertNextCell(NORM_POLYGON, connec) + m1.finishInsertingCells() + + m2 = MEDCouplingDataForTest.buildCircle(0.25, 0.2, 0.4) + # Was looping indefinitly: + m_intersec, resToM1, resToM2 = MEDCouplingUMesh.Intersect2DMeshes(m1, m2, eps) + m_intersec.zipCoords() + coo_tgt = DataArrayDouble([-0.5, -0.5, -0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.03284271247461901, 0.4828427124746191, + -0.014575131106459124, 0.5000000000000001, 0.5, -0.11224989991991996, 0.24271243444677046, 0.5, 0.5, 0.19387505004004, + -0.04799910280454185, -0.06682678787499614, -0.023843325638122054, 0.4915644577163915, 0.5, -0.30612494995996, 0.0, -0.5, + -0.5, 0.0, -0.25728756555322957, 0.5, -0.023843325638122026, 0.49156445771639157, -0.04799910280454181, -0.06682678787499613], 17 ,2) + conn_tgt = [32, 5, 2, 6, 4, 7, 8, 9, 10, 32, 6, 3, 0, 1, 5, 4, 11, 12, 13, 14, 15, 16] + connI_tgt = [0, 9, 22] + res1_tgt = [0, 0] + res2_tgt = [0, -1] + self.assert_(coo_tgt.isEqualWithoutConsideringStr(m_intersec.getCoords(), 1e-12)) + self.assertEqual(conn_tgt, m_intersec.getNodalConnectivity().getValues()) + self.assertEqual(connI_tgt, m_intersec.getNodalConnectivityIndex().getValues()) + self.assertEqual(res1_tgt, resToM1.getValues()) + self.assertEqual(res2_tgt, resToM2.getValues()) + def testDAIBuildUnique1(self): d=DataArrayInt([1,2,2,3,3,3,3,4,5,5,7,7,7,19]) e=d.buildUnique() diff --git a/src/MEDCoupling_Swig/MEDCouplingDataForTest.py b/src/MEDCoupling_Swig/MEDCouplingDataForTest.py index c888a0d74..5ed647cba 100644 --- a/src/MEDCoupling_Swig/MEDCouplingDataForTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingDataForTest.py @@ -683,6 +683,24 @@ class MEDCouplingDataForTest: fff.setGaussLocalizationOnCells([6,7],[-1.,-1.,-1.,-1.,1.,-1.,1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,-1.,1.,1.,1.,1.,1.,1.,-1.,1.],[-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626,-0.577350269189626,-0.577350269189626,0.577350269189626,-0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626,-0.577350269189626,0.577350269189626,0.577350269189626,0.577350269189626],[1.,1.,1.,1.,1.,1.,1.,1.]) return MEDCouplingFieldTemplate(fff) + def buildCircle(self, center_X, center_Y, radius): + from cmath import rect + from math import pi + + c = [rect(radius, i*pi/4.0) for i in range(8)] + coords = [c[-1].real,c[-1].imag, c[3].real,c[3].imag, + c[5].real,c[5].imag, c[1].real,c[1].imag] + connec = range(4) + baseMesh = MEDCouplingUMesh.New("circle", 2) + baseMesh.allocateCells(1) + meshCoords = DataArrayDouble.New(coords, len(coords)/2, 2) + meshCoords += (center_X, center_Y) + baseMesh.setCoords(meshCoords) + + baseMesh.insertNextCell(NORM_QPOLYG, connec) + baseMesh.finishInsertingCells() + return baseMesh + build2DTargetMesh_1=classmethod(build2DTargetMesh_1) build2DSourceMesh_1=classmethod(build2DSourceMesh_1) build3DTargetMesh_1=classmethod(build3DTargetMesh_1) @@ -711,4 +729,8 @@ class MEDCouplingDataForTest: buildFieldOnGauss_2=classmethod(buildFieldOnGauss_2) buildFieldOnGauss_3=classmethod(buildFieldOnGauss_3) buildFieldOnGauss_4=classmethod(buildFieldOnGauss_4) + buildCircle=classmethod(buildCircle) pass + + +