return ret;
}
+int MEDCouplingCMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+ //not implemented yet !
+ return -1;
+}
+
void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
{
throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
MEDCouplingFieldDouble *buildOrthogonalField() const;
+ int getElementContainingPoint(const double *pos, double eps) const;
void rotate(const double *center, const double *vector, double angle);
void translate(const double *vector);
MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField not implemented yet !");
}
+int MEDCouplingExtrudedMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+ //not implemented yet
+ return -1;
+}
+
MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh()
{
if(_mesh2D)
MEDCouplingFieldDouble *getMeasureField(bool) const;
MEDCouplingFieldDouble *getMeasureFieldOnNode(bool) const;
MEDCouplingFieldDouble *buildOrthogonalField() const;
+ int getElementContainingPoint(const double *pos, double eps) const;
static int findCorrespCellByNodalConn(const std::vector<int>& nodalConnec,
const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception);
void rotate(const double *center, const double *vector, double angle);
using namespace ParaMEDMEM;
+const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
+
const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
+MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
+{
+}
+
MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
{
switch(type)
return mesh->getMeasureField(isAbs);
}
+void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+{
+ int id=mesh->getElementContainingPoint(loc,_precision);
+ arr->getTuple(id,res);
+}
+
/*!
* Nothing to do. It's not a bug.
*/
return mesh->getMeasureFieldOnNode(isAbs);
}
+void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
+{
+}
+
void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const
{
int oldNbOfElems=old2New->getNbOfElems();
{
public:
static MEDCouplingFieldDiscretization *New(TypeOfField type);
+ double getPrecision() const { return _precision; }
+ void setPrecision(double val) { _precision=val; }
static TypeOfField getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception);
virtual TypeOfField getEnum() const = 0;
virtual bool isEqual(const MEDCouplingFieldDiscretization *other) const = 0;
virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0;
virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0;
virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0;
+ virtual void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const = 0;
virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0;
virtual void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const = 0;
+ protected:
+ MEDCouplingFieldDiscretization();
+ protected:
+ double _precision;
+ static const double DFLT_PRECISION;
};
class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationP0 : public MEDCouplingFieldDiscretization
void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
+ void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const;
void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const;
MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
public:
void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
+ void getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const;
MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
void renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const;
static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di);
void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception)
{
- _time_discr->checkNoTimePresence();
+ const DataArrayDouble *arr=_time_discr->getArray();
+ _type->getValueOn(arr,_mesh,spaceLoc,res);
}
void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception)
{
- _time_discr->checkTimePresence(time);
+ std::vector< const DataArrayDouble *> arrs=_time_discr->getArraysForTime(time);
+ std::vector<double> res2;
+ for(std::vector< const DataArrayDouble *>::const_iterator iter=arrs.begin();iter!=arrs.end();iter++)
+ {
+ int sz=res2.size();
+ res2.resize(sz+(*iter)->getNumberOfComponents());
+ _type->getValueOn(*iter,_mesh,spaceLoc,&res2[sz]);
+ }
+ _time_discr->getValueForTime(time,res2,res);
}
void MEDCouplingFieldDouble::applyLin(double a, double b, int compoId)
virtual void getBoundingBox(double *bbox) const = 0;
virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0;
virtual MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const = 0;
+ virtual int getElementContainingPoint(const double *pos, double eps) const = 0;
virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const;
virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, const char *func) const;
virtual MEDCouplingFieldDouble *buildOrthogonalField() const = 0;
* 'This' is expected to be of spaceDim==3. Idem for 'center' and 'vect'
*/
void MEDCouplingPointSet::rotate3D(const double *center, const double *vect, double angle)
+{
+ double *coords=_coords->getPointer();
+ int nbNodes=getNumberOfNodes();
+ rotate3DAlg(center,vect,angle,nbNodes,coords);
+}
+
+void MEDCouplingPointSet::rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords)
{
double sina=sin(angle);
double cosa=cos(angle);
std::transform(matrixTmp,matrixTmp+9,matrixTmp,std::bind2nd(std::multiplies<double>(),sina));
std::transform(matrix,matrix+9,matrixTmp,matrix,std::plus<double>());
//rotation matrix computed.
- double *coords=_coords->getPointer();
- int nbNodes=getNumberOfNodes();
double tmp[3];
for(int i=0; i<nbNodes; i++)
{
* 'This' is expected to be of spaceDim==2. Idem for 'center' and 'vect'
*/
void MEDCouplingPointSet::rotate2D(const double *center, double angle)
+{
+ double *coords=_coords->getPointer();
+ int nbNodes=getNumberOfNodes();
+ rotate2DAlg(center,angle,nbNodes,coords);
+}
+
+void MEDCouplingPointSet::rotate2DAlg(const double *center, double angle, int nbNodes, double *coords)
{
double cosa=cos(angle);
double sina=sin(angle);
double matrix[4];
matrix[0]=cosa; matrix[1]=-sina; matrix[2]=sina; matrix[3]=cosa;
- double *coords=_coords->getPointer();
- int nbNodes=getNumberOfNodes();
double tmp[2];
for(int i=0; i<nbNodes; i++)
{
void tryToShareSameCoords(MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception);
static DataArrayDouble *mergeNodesArray(const MEDCouplingPointSet *m1, const MEDCouplingPointSet *m2);
static MEDCouplingPointSet *buildInstanceFromMeshType(MEDCouplingMeshType type);
+ static void rotate2DAlg(const double *center, double angle, int nbNodes, double *coords);
+ static void rotate3DAlg(const double *center, const double *vect, double angle, int nbNodes, double *coords);
virtual MEDCouplingPointSet *buildPartOfMySelf(const int *start, const int *end, bool keepCoords) const = 0;
virtual MEDCouplingPointSet *buildPartOfMySelfNode(const int *start, const int *end, bool fullyIn) const = 0;
virtual void findBoundaryNodes(std::vector<int>& nodes) const = 0;
throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
}
-DataArrayDouble *MEDCouplingNoTimeLabel::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingNoTimeLabel::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
+{
+ throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
+}
+
+void MEDCouplingNoTimeLabel::getValueForTime(double time, const std::vector<double>& vals, double *res) const
{
throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
}
}
}
-DataArrayDouble *MEDCouplingWithTimeStep::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingWithTimeStep::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
{
if(std::fabs(time-_time)<=_time_tolerance)
{
- if(_array)
- _array->incrRef();
- return _array;
+ std::vector< const DataArrayDouble *> ret(1);
+ ret[0]=_array;
+ return ret;
}
else
throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
}
+void MEDCouplingWithTimeStep::getValueForTime(double time, const std::vector<double>& vals, double *res) const
+{
+ std::copy(vals.begin(),vals.end(),res);
+}
+
void MEDCouplingWithTimeStep::getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception)
{
if(std::fabs(time-_time)<=_time_tolerance)
return new MEDCouplingConstOnTimeInterval(*this,deepCpy);
}
-DataArrayDouble *MEDCouplingConstOnTimeInterval::getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception)
+std::vector< const DataArrayDouble *> MEDCouplingConstOnTimeInterval::getArraysForTime(double time) const throw(INTERP_KERNEL::Exception)
{
if(time>_start_time-_time_tolerance && time<_end_time+_time_tolerance)
{
- if(_array)
- _array->incrRef();
- return _array;
+ std::vector< const DataArrayDouble *> ret(1);
+ ret[0]=_array;
+ return ret;
}
else
throw INTERP_KERNEL::Exception(EXCEPTION_MSG);
}
+void MEDCouplingConstOnTimeInterval::getValueForTime(double time, const std::vector<double>& vals, double *res) const
+{
+ std::copy(vals.begin(),vals.end(),res);
+}
+
bool MEDCouplingConstOnTimeInterval::areCompatible(const MEDCouplingTimeDiscretization *other) const
{
if(!MEDCouplingTimeDiscretization::areCompatible(other))
virtual void setArrays(const std::vector<DataArrayDouble *>& arrays, TimeLabel *owner) throw(INTERP_KERNEL::Exception);
DataArrayDouble *getArray() const { return _array; }
virtual DataArrayDouble *getEndArray() const { return _array; }
- //! Warning contrary to getArray method this method returns an object to deal with.
- virtual DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception) = 0;
+ virtual void getValueForTime(double time, const std::vector<double>& vals, double *res) const = 0;
virtual void getArrays(std::vector<DataArrayDouble *>& arrays) const;
virtual bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
virtual bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
void checkNoTimePresence() const throw(INTERP_KERNEL::Exception) { }
void checkTimePresence(double time) const throw(INTERP_KERNEL::Exception);
- DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+ std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+ void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
bool isBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
bool isStrictlyBefore(const MEDCouplingTimeDiscretization *other) const throw(INTERP_KERNEL::Exception);
double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception);
double getStartTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
double getEndTime(int& dt, int& it) const throw(INTERP_KERNEL::Exception) { dt=_dt; it=_it; return _time; }
DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+ std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+ void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
public:
MEDCouplingTimeDiscretization *performCpy(bool deepCpy) const;
bool areCompatible(const MEDCouplingTimeDiscretization *other) const;
bool isEqual(const MEDCouplingTimeDiscretization *other, double prec) const;
- DataArrayDouble *getArrayOnTime(double time) const throw(INTERP_KERNEL::Exception);
+ std::vector< const DataArrayDouble *> getArraysForTime(double time) const throw(INTERP_KERNEL::Exception);
+ void getValueForTime(double time, const std::vector<double>& vals, double *res) const;
void getValueOnTime(int eltId, double time, double *value) const throw(INTERP_KERNEL::Exception);
void getValueOnDiscTime(int eltId, int dt, int it, double *value) const throw(INTERP_KERNEL::Exception);
TypeOfTimeDiscretization getEnum() const { return DISCRETIZATION; }
#include "CellModel.hxx"
#include "VolSurfUser.txx"
#include "InterpolationUtils.hxx"
+#include "PointLocatorAlgos.txx"
+#include "BBTree.txx"
#include <sstream>
#include <limits>
return ret;
}
+/*{
+ const double *coordsPtr=_coords->getConstPointer();
+ BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbNodes,-prec);
+ double bb[2*SPACEDIM];
+ double prec2=prec*prec;
+ for(int i=0;i<nbNodes;i++)
+ {
+ if(std::find(c.begin(),c.end(),i)!=c.end())
+ continue;
+ for(int j=0;j<SPACEDIM;j++)
+ {
+ bb[2*j]=coordsPtr[SPACEDIM*i+j];
+ bb[2*j+1]=coordsPtr[SPACEDIM*i+j];
+ }
+ std::vector<int> intersectingElems;
+ myTree.getIntersectingElems(bb,intersectingElems);
+ if(intersectingElems.size()>1)
+ {
+ std::vector<int> commonNodes;
+ for(std::vector<int>::const_iterator it=intersectingElems.begin();it!=intersectingElems.end();it++)
+ if(*it!=i)
+ if(INTERP_KERNEL::distance2<SPACEDIM>(coordsPtr+SPACEDIM*i,coordsPtr+SPACEDIM*(*it))<prec2)
+ commonNodes.push_back(*it);
+ if(!commonNodes.empty())
+ {
+ cI.push_back(cI.back()+commonNodes.size()+1);
+ c.push_back(i);
+ c.insert(c.end(),commonNodes.begin(),commonNodes.end());
+ }
+ }
+ }*/
+
+int MEDCouplingUMesh::getElementContainingPoint(const double *pos, double eps) const
+{
+ std::vector<int> elts;
+ getElementsContainingPoint(pos,eps,elts);
+ if(elts.empty())
+ return -1;
+ return elts.front();
+}
+
+void MEDCouplingUMesh::getElementsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const
+{
+ std::vector<int> eltsIndex;
+ getElementsContainingPoints(pos,1,eps,elts,eltsIndex);
+}
+
+namespace ParaMEDMEM
+{
+ template<const int SPACEDIMM>
+ class DummyClsMCUG
+ {
+ public:
+ static const int MY_SPACEDIM=SPACEDIMM;
+ static const int MY_MESHDIM=8;
+ typedef int MyConnType;
+ static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
+ };
+}
+
+template<int SPACEDIM>
+void MEDCouplingUMesh::getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
+ double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const
+{
+ std::vector<double> bbox;
+ eltsIndex.resize(nbOfPoints+1);
+ eltsIndex[0]=0;
+ elts.clear();
+ getBoundingBoxForBBTree(bbox);
+ int nbOfCells=getNumberOfCells();
+ const int *conn=_nodal_connec->getConstPointer();
+ const int *connI=_nodal_connec_index->getConstPointer();
+ double bb[2*SPACEDIM];
+ for(int j=0;j<SPACEDIM;j++)
+ { bb[2*j]=std::numeric_limits<double>::max(); bb[2*j+1]=-std::numeric_limits<double>::max(); }
+ BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbOfCells,-eps);
+ for(int i=0;i<nbOfPoints;i++)
+ {
+ eltsIndex[i+1]=eltsIndex[i];
+ for(int j=0;j<SPACEDIM;j++)
+ {
+ bb[2*j]=std::min(bb[2*j],pos[SPACEDIM*i+j]);
+ bb[2*j+1]=std::max(bb[2*j+1],pos[SPACEDIM*i+j]);
+ }
+ std::vector<int> candidates;
+ myTree.getIntersectingElems(bb,candidates);
+ for(std::vector<int>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
+ {
+ int sz=connI[(*iter)+1]-connI[*iter]-1;
+ if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,
+ (INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]],
+ coords,conn+connI[*iter]+1,sz,eps))
+ {
+ eltsIndex[i+1]++;
+ elts.push_back(*iter);
+ }
+ }
+ }
+}
+
+void MEDCouplingUMesh::getElementsContainingPoints(const double *pos, int nbOfPoints, double eps,
+ std::vector<int>& elts, std::vector<int>& eltsIndex) const
+{
+ int spaceDim=getSpaceDimension();
+ int mDim=getMeshDimension();
+ if(spaceDim==3)
+ {
+ if(mDim==3)
+ {
+ const double *coords=_coords->getConstPointer();
+ getElementsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
+ }
+ /*else if(mDim==2)
+ {
+
+ }*/
+ else
+ throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
+ }
+ else if(spaceDim==2)
+ {
+ if(mDim==2)
+ {
+ const double *coords=_coords->getConstPointer();
+ getElementsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
+ }
+ else
+ throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !");
+ }
+
+
+
+}
+
/*!
* This method is only available for a mesh with meshDim==2 and spaceDim==2||spaceDim==3.
* This method returns a vector 'cells' where all detected butterfly cells have been added to cells.
}
}
+/*!
+ * This method aggregate the bbox of each cell and put it into bbox parameter.
+ * @param bbox out parameter of size 2*spacedim*nbOfcells.
+ */
+void MEDCouplingUMesh::getBoundingBoxForBBTree(std::vector<double>& bbox) const
+{
+ int spaceDim=getSpaceDimension();
+ int nbOfCells=getNumberOfCells();
+ bbox.resize(2*nbOfCells*spaceDim);
+ for(int i=0;i<nbOfCells*spaceDim;i++)
+ {
+ bbox[2*i]=std::numeric_limits<double>::max();
+ bbox[2*i+1]=-std::numeric_limits<double>::max();
+ }
+ const double *coordsPtr=_coords->getConstPointer();
+ const int *conn=_nodal_connec->getConstPointer();
+ const int *connI=_nodal_connec_index->getConstPointer();
+ for(int i=0;i<nbOfCells;i++)
+ {
+ int offset=connI[i]+1;
+ int nbOfNodesForCell=connI[i+1]-offset;
+ for(int j=0;j<nbOfNodesForCell;j++)
+ {
+ int nodeId=conn[offset+j];
+ if(nodeId>=0)
+ for(int k=0;k<spaceDim;k++)
+ {
+ bbox[2*spaceDim*i+2*k]=std::min(bbox[2*spaceDim*i+2*k],coordsPtr[spaceDim*nodeId+k]);
+ bbox[2*spaceDim*i+2*k+1]=std::max(bbox[2*spaceDim*i+2*k+1],coordsPtr[spaceDim*nodeId+k]);
+ }
+ }
+ }
+}
+
namespace ParaMEDMEMImpl
{
class ConnReader
MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
MEDCouplingFieldDouble *buildOrthogonalField() const;
+ int getElementContainingPoint(const double *pos, double eps) const;
+ void getElementsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const;
+ void getElementsContainingPoints(const double *pos, int nbOfPoints, double eps,
+ std::vector<int>& elts, std::vector<int>& eltsIndex) const;
void checkButterflyCells(std::vector<int>& cells) const;
+ void getBoundingBoxForBBTree(std::vector<double>& bbox) const;
bool checkConsecutiveCellTypes() const;
MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
DataArrayDouble *getBarycenterAndOwner() const;
void checkFullyDefined() const throw(INTERP_KERNEL::Exception);
//tools
MEDCouplingUMesh *buildPartOfMySelfKeepCoords(const int *start, const int *end) const;
+ template<int SPACEDIM>
+ void getElementsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
+ double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const;
private:
//! this iterator stores current position in _nodal_connec array.
mutable int _iterator;
//not implemented yet.
return 0;
}
+
+int MEDCouplingUMeshDesc::getElementContainingPoint(const double *pos, double eps) const
+{
+ //not implemented yet.
+ return -1;
+}
void renumberNodes(const int *newNodeNumbers, int newNbOfNodes);
MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
+ int getElementContainingPoint(const double *pos, double eps) const;
MEDCouplingFieldDouble *buildOrthogonalField() const;
DataArrayInt *zipCoordsTraducer();
MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
targetMesh->decrRef();
}
+void MEDCouplingBasicsTest::testGetElementsContainingPoint()
+{
+ MEDCouplingUMesh *targetMesh=build2DTargetMesh_1();
+ double pos[12]={0.,0.,0.4,0.4,0.,0.4,0.1,0.1,0.25,0.,0.65,0.};
+ std::vector<int> t1,t2;
+ //2D basic
+ targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2);
+ CPPUNIT_ASSERT_EQUAL(6,(int)t1.size());
+ CPPUNIT_ASSERT_EQUAL(7,(int)t2.size());
+ const int expectedValues1[6]={0,4,3,0,1,2};
+ const int expectedValues2[7]={0,1,2,3,4,5,6};
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1));
+ CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2));
+ //2D with no help of bounding box.
+ double center[2]={0.2,0.2};
+ MEDCouplingPointSet::rotate2DAlg(center,0.78539816339744830962,6,pos);
+ targetMesh->rotate(center,0,0.78539816339744830962);
+ t1.clear(); t2.clear();
+ targetMesh->getElementsContainingPoints(pos,6,1e-12,t1,t2);
+ CPPUNIT_ASSERT_EQUAL(6,(int)t1.size());
+ CPPUNIT_ASSERT_EQUAL(7,(int)t2.size());
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues1));
+ CPPUNIT_ASSERT(std::equal(t2.begin(),t2.end(),expectedValues2));
+ //2D outside
+ const double pos1bis[2]={-0.3303300858899107,-0.11819805153394641};
+ CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos1bis,1e-12));
+ targetMesh->decrRef();
+ //test limits 2D
+ targetMesh=build2DTargetMesh_1();
+ const double pos2[2]={0.2,-0.05};
+ t1.clear();
+ targetMesh->getElementsContainingPoint(pos2,1e-12,t1);
+ CPPUNIT_ASSERT_EQUAL(2,(int)t1.size());
+ const int expectedValues3[2]={0,1};
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues3));
+ const double pos3[2]={0.2,0.2};
+ t1.clear();
+ targetMesh->getElementsContainingPoint(pos3,1e-12,t1);
+ CPPUNIT_ASSERT_EQUAL(5,(int)t1.size());
+ const int expectedValues4[5]={0,1,2,3,4};
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues4));
+ CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos3,1e-12));
+ targetMesh->decrRef();
+ //3D
+ targetMesh=build3DTargetMesh_1();
+ const double pos4[3]={25.,25.,25.};
+ CPPUNIT_ASSERT_EQUAL(0,targetMesh->getElementContainingPoint(pos4,1e-12));
+ const double pos5[3]={50.,50.,50.};
+ t1.clear();
+ targetMesh->getElementsContainingPoint(pos5,1e-12,t1);
+ CPPUNIT_ASSERT_EQUAL(8,(int)t1.size());
+ const int expectedValues5[8]={0,1,2,3,4,5,6,7};
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues5));
+ const double pos6[3]={0., 50., 0.};
+ t1.clear();
+ targetMesh->getElementsContainingPoint(pos6,1e-12,t1);
+ CPPUNIT_ASSERT_EQUAL(2,(int)t1.size());
+ const int expectedValues6[2]={0,2};
+ CPPUNIT_ASSERT(std::equal(t1.begin(),t1.end(),expectedValues6));
+ //3D outside
+ const double pos7[3]={-1.0,-1.0,0.};
+ CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos7,1e-12));
+ //3D outside 2
+ const double center2[3]={0.,0.,0.};
+ const double vec2[3]={0.,-1.,0.};
+ targetMesh->rotate(center2,vec2,0.78539816339744830962);
+ const double pos8[3]={-25,25.,12.};
+ CPPUNIT_ASSERT_EQUAL(-1,targetMesh->getElementContainingPoint(pos8,1e-12));
+ //
+ targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::testGetValueOn1()
+{
+ MEDCouplingUMesh *targetMesh=build2DTargetMesh_1();
+ MEDCouplingFieldDouble *fieldOnCells=MEDCouplingFieldDouble::New(ON_CELLS);
+ int nbOfCells=targetMesh->getNumberOfCells();
+ fieldOnCells->setMesh(targetMesh);
+ DataArrayDouble *array=DataArrayDouble::New();
+ array->alloc(nbOfCells,2);
+ fieldOnCells->setArray(array);
+ double *tmp=array->getPointer();
+ for(int i=0;i<nbOfCells;i++)
+ { tmp[2*i]=7.+(double)i; tmp[2*i+1]=17.+(double)i; }
+ array->decrRef();
+ //
+ const double pos1[2]={0.25,0.};
+ double res[2];
+ fieldOnCells->getValueOn(pos1,res);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(8.,res[0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(18.,res[1],1e-12);
+ //
+ fieldOnCells->decrRef();
+ targetMesh->decrRef();
+}
+
void MEDCouplingBasicsTest::test2DInterpP0P0_1()
{
MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1();
CPPUNIT_TEST( testMergeNodesOnField );
CPPUNIT_TEST( testCheckConsecutiveCellTypes );
CPPUNIT_TEST( testBuildOrthogonalField );
+ CPPUNIT_TEST( testGetElementsContainingPoint );
+ CPPUNIT_TEST( testGetValueOn1 );
CPPUNIT_TEST( test2DInterpP0P0_1 );
CPPUNIT_TEST( test2DInterpP0P0PL_1 );
CPPUNIT_TEST( test2DInterpP0P0PL_2 );
void testMergeNodesOnField();
void testCheckConsecutiveCellTypes();
void testBuildOrthogonalField();
+ void testGetElementsContainingPoint();
+ void testGetValueOn1();
void test2DInterpP0P0_1();
void test2DInterpP0P0PL_1();
void test2DInterpP0P0PL_2();
