QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
{
- QuadraticPolygon *ret=new QuadraticPolygon;
+ QuadraticPolygon *ret(new QuadraticPolygon);
std::size_t size=nodes.size();
for(std::size_t i=0;i<size;i++)
{
QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
{
- QuadraticPolygon *ret=new QuadraticPolygon;
+ QuadraticPolygon *ret(new QuadraticPolygon);
std::size_t size=nodes.size();
for(std::size_t i=0;i<size/2;i++)
+
{
EdgeLin *e1,*e2;
e1=new EdgeLin(nodes[i],nodes[i+size/2]);
return ret;
}
+Edge *QuadraticPolygon::BuildLinearEdge(std::vector<Node *>& nodes)
+{
+ if(nodes.size()!=2)
+ throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildLinearEdge : input vector is expected to be of size 2 !");
+ Edge *ret(new EdgeLin(nodes[0],nodes[1]));
+ nodes[0]->decrRef(); nodes[1]->decrRef();
+ return ret;
+}
+
+Edge *QuadraticPolygon::BuildArcCircleEdge(std::vector<Node *>& nodes)
+{
+ if(nodes.size()!=3)
+ throw INTERP_KERNEL::Exception("QuadraticPolygon::BuildArcCircleEdge : input vector is expected to be of size 3 !");
+ EdgeLin *e1(new EdgeLin(nodes[0],nodes[2])),*e2(new EdgeLin(nodes[2],nodes[1]));
+ SegSegIntersector inters(*e1,*e2);
+ bool colinearity=inters.areColinears();
+ delete e1; delete e2;
+ Edge *ret(0);
+ if(colinearity)
+ ret=new EdgeLin(nodes[0],nodes[1]);
+ else
+ ret=new EdgeArcCircle(nodes[0],nodes[2],nodes[1]);
+ nodes[0]->decrRef(); nodes[1]->decrRef(); nodes[2]->decrRef();
+ return ret;
+}
+
void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
{
std::ofstream file(fileName);
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 Edge *BuildLinearEdge(std::vector<Node *>& nodes);
+ INTERPKERNEL_EXPORT static Edge *BuildArcCircleEdge(std::vector<Node *>& nodes);
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((mDim==3 && sDim==3) || (mDim==2 && sDim==3) || (mDim==1 && sDim==1) || ( mDim==1 && sDim==3)) // Compute refined boundary box for quadratic elements only in 2D.
return getBoundingBoxForBBTreeFast();
- else
+ if((mDim==2 && sDim==2) || (mDim==1 && sDim==2))
{
bool presenceOfQuadratic(false);
for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=_types.begin();it!=_types.end();it++)
if(cm.isQuadratic())
presenceOfQuadratic=true;
}
- if(presenceOfQuadratic)
+ if(!presenceOfQuadratic)
+ return getBoundingBoxForBBTreeFast();
+ if(mDim==2 && sDim==2)
return getBoundingBoxForBBTree2DQuadratic(arcDetEps);
else
- return getBoundingBoxForBBTreeFast();
+ return getBoundingBoxForBBTree1DQuadratic(arcDetEps);
}
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree : Managed dimensions are (mDim=1,sDim=1), (mDim=1,sDim=2), (mDim=1,sDim=3), (mDim=2,sDim=2), (mDim=2,sDim=3) and (mDim=3,sDim=3) !");
}
/*!
}
/*!
- * 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.
* \throw If \a this is not fully defined.
* \throw If \a this is not a mesh with meshDimension equal to 2.
* \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic
*/
DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
{
checkFullyDefined();
int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
- if(mDim!=2 || spaceDim!=2)
+ if(spaceDim!=2 || 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());
{
const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
int sz(connI[1]-connI[0]-1);
- INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1e-12;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
std::vector<INTERP_KERNEL::Node *> nodes(sz);
INTERP_KERNEL::QuadraticPolygon *pol(0);
for(int j=0;j<sz;j++)
return ret.retn();
}
+/*!
+ * This method aggregates the bbox of each 1D 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.
+ * \throw If \a this is not fully defined.
+ * \throw If \a this is not a mesh with meshDimension equal to 1.
+ * \throw If \a this is not a mesh with spaceDimension equal to 2.
+ * \sa MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic
+ */
+DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arcDetEps) const
+{
+ checkFullyDefined();
+ int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells()),nbOfNodes(getNumberOfNodes());
+ if(spaceDim!=2 || spaceDim!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic : This method should be applied on mesh with mesh dimension equal to 1 and space dimension also equal to 2!");
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
+ double *bbox(ret->getPointer());
+ const double *coords(_coords->getConstPointer());
+ const int *conn(_nodal_connec->getConstPointer()),*connI(_nodal_connec_index->getConstPointer());
+ for(int i=0;i<nbOfCells;i++,bbox+=4,connI++)
+ {
+ const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
+ int sz(connI[1]-connI[0]-1);
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
+ std::vector<INTERP_KERNEL::Node *> nodes(sz);
+ INTERP_KERNEL::Edge *edge(0);
+ for(int j=0;j<sz;j++)
+ {
+ int nodeId(conn[*connI+1+j]);
+ nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*2],coords[nodeId*2+1]);
+ }
+ if(!cm.isQuadratic())
+ edge=INTERP_KERNEL::QuadraticPolygon::BuildLinearEdge(nodes);
+ else
+ edge=INTERP_KERNEL::QuadraticPolygon::BuildArcCircleEdge(nodes);
+ const INTERP_KERNEL::Bounds& b(edge->getBounds());
+ bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); edge->decrRef();
+ }
+ return ret.retn();
+}
+
/// @cond INTERNAL
namespace ParaMEDMEMImpl
MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree(double arcDetEps=1e-12) const;
MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTreeFast() const;
MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const;
+ MEDCOUPLING_EXPORT DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const;
MEDCOUPLING_EXPORT MEDCouplingUMesh *buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy);
MEDCOUPLING_EXPORT bool isFullyQuadratic() const;
MEDCOUPLING_EXPORT bool isPresenceOfQuadratic() const;
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()
%newobject ParaMEDMEM::MEDCouplingUMesh::buildUnionOf3DMesh;
%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTreeFast;
%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic;
+%newobject ParaMEDMEM::MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic;
%newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__;
%newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__;
%newobject ParaMEDMEM::MEDCoupling1GTUMesh::New;
DataArrayInt *buildUnionOf3DMesh() const throw(INTERP_KERNEL::Exception);
DataArrayDouble *getBoundingBoxForBBTreeFast() const throw(INTERP_KERNEL::Exception);
DataArrayDouble *getBoundingBoxForBBTree2DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
+ DataArrayDouble *getBoundingBoxForBBTree1DQuadratic(double arcDetEps=1e-12) const throw(INTERP_KERNEL::Exception);
static MEDCouplingUMesh *Build0DMeshFromCoords(DataArrayDouble *da) throw(INTERP_KERNEL::Exception);
static MEDCouplingUMesh *MergeUMeshes(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception);
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
+
+
+
