#include "InterpKernelGeo2DQuadraticPolygon.hxx"
#include "InterpKernelGeo2DNode.hxx"
#include "DirectedBoundingBox.hxx"
+#include "InterpKernelAutoPtr.hxx"
#include <cmath>
#include <limits>
void MEDCouplingPointSet::findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
{
if(getSpaceDimension()!=3)
- throw INTERP_KERNEL::Exception("Invalid spacedim to be applied on this ! Must be equal to 3 !");
+ throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnPlane : Invalid spacedim to be applied on this ! Must be equal to 3 !");
int nbOfNodes=getNumberOfNodes();
double a=vec[0],b=vec[1],c=vec[2],d=-pt[0]*vec[0]-pt[1]*vec[1]-pt[2]*vec[2];
double deno=sqrt(a*a+b*b+c*c);
}
}
+/*!
+ * This method is expecting to be called for meshes so that getSpaceDimension() returns 2 or 3.
+ * This method returns in 'nodes' output all the nodes that are at a distance lower than epsilon from a line
+ * defined by the point 'pt' and the vector 'vec'. 'pt' and 'vec' are expected to have a dimension equal to space dimension of 'this'
+ * @param pt points to an array of size this->getSpaceDimension and represents a point that owns to plane.
+ * @param vec points to an array of size this->getSpaceDimension and represents the direction vector of the line. The norm of the vector is not compulsory equal to 1.
+ * But norm must be greater than 10*abs(eps)
+ * @param eps is the maximal distance around the plane where node in this->_coords will be picked.
+ * @param nodes is the output of the method. The vector is not compulsory empty before call. The nodes that fulfills the condition will be added at the end of the nodes.
+ */
+void MEDCouplingPointSet::findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception)
+{
+ int spaceDim=getSpaceDimension();
+ if(spaceDim!=2 && spaceDim!=3)
+ throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid spacedim to be applied on this ! Must be equal to 2 or 3 !");
+ int nbOfNodes=getNumberOfNodes();
+ int den=0.;
+ for(int i=0;i<spaceDim;i++)
+ den+=vec[i]*vec[i];
+ double deno=sqrt(den);
+ if(deno<10.*eps)
+ throw INTERP_KERNEL::Exception("MEDCouplingPointSet::findNodesOnLine : Invalid given direction vector ! Norm is too small !");
+ INTERP_KERNEL::AutoPtr<double> vecn=new double[spaceDim];
+ for(int i=0;i<spaceDim;i++)
+ vecn[i]=vec[i]/deno;
+ const double *work=_coords->getConstPointer();
+ if(spaceDim==2)
+ {
+ for(int i=0;i<nbOfNodes;i++)
+ {
+ if(std::abs(vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]))<eps)
+ nodes.push_back(i);
+ work+=2;
+ }
+ }
+ else
+ {
+ for(int i=0;i<nbOfNodes;i++)
+ {
+ double a=vecn[0]*(work[1]-pt[1])-vecn[1]*(work[0]-pt[0]);
+ double b=vecn[1]*(work[2]-pt[2])-vecn[2]*(work[1]-pt[1]);
+ double c=vecn[2]*(work[0]-pt[0])-vecn[0]*(work[2]-pt[2]);
+ if(std::sqrt(a*a+b*b+c*c)<eps)
+ nodes.push_back(i);
+ work+=3;
+ }
+ }
+}
+
/*!
* merge _coords arrays of m1 and m2 and returns the union. The returned instance is newly created with ref count == 1.
*/
void tryToShareSameCoords(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception);
virtual void tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception) = 0;
void findNodesOnPlane(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception);
+ void findNodesOnLine(const double *pt, const double *vec, double eps, std::vector<int>& nodes) const throw(INTERP_KERNEL::Exception);
static DataArrayDouble *MergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2) throw(INTERP_KERNEL::Exception);
static DataArrayDouble *MergeNodesArray(const std::vector<const MEDCouplingPointSet *>& ms) throw(INTERP_KERNEL::Exception);
static MEDCouplingPointSet *BuildInstanceFromMeshType(MEDCouplingMeshType type);
std::map<int,INTERP_KERNEL::Node *> mappRev;
INTERP_KERNEL::QuadraticPolygon pol1;
MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,b1[i],desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */pol1,mapp,mappRev);
+ std::vector<int> crTmp,crITmp;
+ crITmp.push_back(crI.back());
for(std::vector<int>::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++)
{
INTERP_KERNEL::QuadraticPolygon pol2;
MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,b2[*it2],desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */pol2,mapp,mappRev);
pol1.buildPartitionsAbs(pol2,mapp,i,*it2,cr,crI,cNb1,cNb2);
}
+ if(!crTmp.empty())
+ {
+ cr.insert(cr.end(),crTmp.begin(),crTmp.end());
+ crI.insert(crI.end(),crITmp.begin()+1,crITmp.end());
+ }
for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
(*it).second->decrRef();
}
}
}
+/*!
+ * Given a 2D mesh conn by (conn2D,connI2D) it returns a single polygon
+ */
+void MEDCouplingUMesh::BuildUnionOf2DMesh(const std::vector<int>& conn2D, const std::vector<int>& connI2D, std::vector<int>& polyUnion)
+{
+}
+
MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)),
_own_cell(true),_cell_id(-1),_nb_cell(0)
{
const MEDCouplingUMesh *m2, const std::vector<bool>& b2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
const std::vector<double>& addCoords,
std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2);
+ static void BuildUnionOf2DMesh(const std::vector<int>& conn2D, const std::vector<int>& connI2D, std::vector<int>& polyUnion);
/// @endcond
private:
//! this iterator stores current position in _nodal_connec array.
coordY->decrRef();
m1c->decrRef();
}
+
+void MEDCouplingBasicsTest4::testFindNodesOnLine1()
+{
+ MEDCouplingUMesh *mesh=build2DTargetMesh_1();
+ const double pt[2]={-0.3,-0.3};
+ const double pt2[3]={0.,0.,0.};
+ const double pt3[3]={-0.3,0.,0.};
+ const double vec[2]={0.,1.};
+ const double vec2[3]={1.,0.,0.};
+ const double vec3[3]={0.,1.,1.};
+ const int expected1[3]={0,3,6};
+ std::vector<int> res;
+ mesh->findNodesOnLine(pt,vec,1e-12,res);
+ CPPUNIT_ASSERT_EQUAL(3,(int)res.size());
+ CPPUNIT_ASSERT(std::equal(expected1,expected1+3,res.begin()));
+ res.clear();
+ //
+ mesh->changeSpaceDimension(3);
+ mesh->rotate(pt2,vec2,M_PI/4.);
+ mesh->findNodesOnLine(pt3,vec3,1e-12,res);
+ CPPUNIT_ASSERT_EQUAL(3,(int)res.size());
+ CPPUNIT_ASSERT(std::equal(expected1,expected1+3,res.begin()));
+ //
+ mesh->decrRef();
+}
CPPUNIT_TEST( testGetNodeIdsInUse1 );
CPPUNIT_TEST( testBuildDescendingConnec2 );
CPPUNIT_TEST( testIntersect2DMeshesTmp1 );
+ CPPUNIT_TEST( testFindNodesOnLine1 );
CPPUNIT_TEST_SUITE_END();
public:
void testDescriptionInMeshTimeUnit1();
void testGetNodeIdsInUse1();
void testBuildDescendingConnec2();
void testIntersect2DMeshesTmp1();
+ void testFindNodesOnLine1();
};
}
self->renumberNodes2(da2->getConstPointer(),newNbOfNodes);
}
}
+ PyObject *findNodesOnLine(PyObject *pt, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception)
+ {
+ std::vector<int> nodes;
+ int spaceDim=self->getSpaceDimension();
+ int sz1,sz2;
+ double *p=convertPyToNewDblArr2(pt,&sz1);
+ double *v=convertPyToNewDblArr2(vec,&sz2);
+ if(sz1!=spaceDim || sz2!=spaceDim)
+ throw INTERP_KERNEL::Exception("Mismatch of spaceDimension and the length of the input array point and vector !");
+ self->findNodesOnLine(p,v,eps,nodes);
+ delete [] v;
+ delete [] p;
+ return convertIntArrToPyList2(nodes);
+ }
PyObject *findNodesOnPlane(PyObject *pt, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception)
{
std::vector<int> nodes;
self.assertAlmostEqual(expected5[i],m3.getCoords().getIJ(0,i),12);
pass
pass
+
+ def testFindNodesOnLine1(self):
+ mesh=MEDCouplingDataForTest.build2DTargetMesh_1();
+ pt=[-0.3,-0.3]
+ pt2=[0.,0.,0.]
+ pt3=[-0.3,0.,0.]
+ vec=[0.,1.]
+ vec2=[1.,0.,0.]
+ vec3=[0.,1.,1.]
+ expected1=[0,3,6]
+ res=mesh.findNodesOnLine(pt,vec,1e-12);
+ self.assertEqual(3,len(res));
+ self.assertEqual(expected1,res);
+ #
+ mesh.changeSpaceDimension(3);
+ mesh.rotate(pt2,vec2,pi/4.);
+ res=mesh.findNodesOnLine(pt3,vec3,1e-12);
+ self.assertEqual(3,len(res));
+ self.assertEqual(expected1,res);
+ pass
def setUp(self):
pass
