was not properly computed. Was leading to spurious bugs in Intersect2DMeshes().
{
}
-QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
+QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes, bool isClosed)
{
QuadraticPolygon *ret=new QuadraticPolygon;
std::size_t size=nodes.size();
- for(std::size_t i=0;i<size;i++)
+ for(std::size_t i=0;i< (size - (isClosed?0:1));i++)
{
ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
nodes[i]->decrRef();
}
+ if(!isClosed)
+ nodes[size-1]->decrRef();
return ret;
}
-QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
+QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes, bool isClosed)
{
QuadraticPolygon *ret=new QuadraticPolygon;
std::size_t size=nodes.size();
- for(std::size_t i=0;i<size/2;i++)
+ std::size_t quad_offset = isClosed ? (size/2) : (size/2+1);
+ for(std::size_t i = 0; i < size/2; i++)
{
EdgeLin *e1,*e2;
- e1=new EdgeLin(nodes[i],nodes[i+size/2]);
- e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]);
+ e1=new EdgeLin(nodes[i],nodes[i+quad_offset]);
+ e2=new EdgeLin(nodes[i+quad_offset],nodes[(i+1)%quad_offset]);
SegSegIntersector inters(*e1,*e2);
bool colinearity=inters.areColinears();
delete e1; delete e2;
if(colinearity)
- ret->pushBack(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;
}
INTERPKERNEL_EXPORT QuadraticPolygon() { }
INTERPKERNEL_EXPORT QuadraticPolygon(const QuadraticPolygon& other):ComposedEdge(other) { }
INTERPKERNEL_EXPORT QuadraticPolygon(const char *fileName);
- INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector<Node *>& nodes);
- INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector<Node *>& nodes);
+ INTERPKERNEL_EXPORT static QuadraticPolygon *BuildLinearPolygon(std::vector<Node *>& nodes, bool isClosed=true);
+ INTERPKERNEL_EXPORT static QuadraticPolygon *BuildArcCirclePolygon(std::vector<Node *>& nodes, bool isClosed=true);
INTERPKERNEL_EXPORT static void BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName);
INTERPKERNEL_EXPORT ~QuadraticPolygon();
INTERPKERNEL_EXPORT void closeMe() 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
{
}
/*!
- * 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.
{
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<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
double *bbox(ret->getPointer());
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();
}
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};
CPPUNIT_TEST( testDAIPartitionByDifferentValues1 );
CPPUNIT_TEST( testDAICheckMonotonic1 );
CPPUNIT_TEST( testIntersect2DMeshesTmp6 );
+ CPPUNIT_TEST( testIntersect2DMeshesTmp7 );
CPPUNIT_TEST( testDAIBuildSubstractionOptimized1 );
CPPUNIT_TEST( testDAIIsStrictlyMonotonic1 );
CPPUNIT_TEST( testSimplexize3 );
void testDAIPartitionByDifferentValues1();
void testDAICheckMonotonic1();
void testIntersect2DMeshesTmp6();
+ void testIntersect2DMeshesTmp7();
void testDAIBuildSubstractionOptimized1();
void testDAIIsStrictlyMonotonic1();
void testSimplexize3();
# 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()
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)
buildFieldOnGauss_2=classmethod(buildFieldOnGauss_2)
buildFieldOnGauss_3=classmethod(buildFieldOnGauss_3)
buildFieldOnGauss_4=classmethod(buildFieldOnGauss_4)
+ buildCircle=classmethod(buildCircle)
pass
+
+
+
