From f0e9f1e8bd80fe9c81987138e5b53c0528e2f151 Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 17 Feb 2012 16:14:40 +0000 Subject: [PATCH] getCellIdsCrossingPlane implementation. --- src/MEDCoupling/MEDCouplingPointSet.hxx | 2 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 66 +++++++++++++++++++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 4 +- src/MEDCoupling/MEDCouplingUMeshDesc.cxx | 2 +- src/MEDCoupling/MEDCouplingUMeshDesc.hxx | 2 +- src/MEDCoupling_Swig/MEDCoupling.i | 32 +++++++++++- 6 files changed, 102 insertions(+), 6 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingPointSet.hxx b/src/MEDCoupling/MEDCouplingPointSet.hxx index e84d46fa5..0cb17a662 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.hxx +++ b/src/MEDCoupling/MEDCouplingPointSet.hxx @@ -102,7 +102,7 @@ namespace ParaMEDMEM void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; void unserialization(const std::vector& tinyInfoD, const std::vector& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); - virtual void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) = 0; + virtual void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) const = 0; virtual void getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems) = 0; virtual DataArrayInt *zipCoordsTraducer() = 0; protected: diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index cca5a8fe1..0efbe6192 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -1595,7 +1595,7 @@ void MEDCouplingUMesh::renumberCells(const int *old2NewBg, bool check) throw(INT * Warning 'elems' is incremented during the call so if elems is not empty before call returned elements will be * added in 'elems' parameter. */ -void MEDCouplingUMesh::getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) +void MEDCouplingUMesh::getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) const { if(getMeshDimension()==-1) { @@ -2399,6 +2399,70 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildDirectionVectorField() const return ret; } +/*! + * This method expects that 'this' is fully defined and has a spaceDim==3 and a meshDim==3. If it is not the case an exception will be thrown. + * This method returns 2 objects : + * - a newly created mesh instance containing the result of the slice lying on the same coords than 'this' and with a meshdim == 2 + * - a newly created dataarray having number of tuples equal to the number of cells in returned mesh that tells for each 2D cell in returned + * mesh the 3D cell id is 'this' it comes from. + * @param origin is the origin of the plane. It should be an array of length 3. + * @param vec is the direction vector of the plane. It should be an array of length 3. Norm of 'vec' should be > 1e-6. + * @param eps is the precision. It is used by called method MEDCouplingUMesh::getCellIdsCrossingPlane for the first 3D cell selection. 'eps' is + * also used to state if the plane shares 2 3D cells. In this case an exception will be thrown. + */ +MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const double *vec, double eps, DataArrayInt *&cellIds) const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=3 || getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!"); + MEDCouplingAutoRefCountObjectPtr candidates=getCellIdsCrossingPlane(origin,vec,eps); + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : not implemented yet !"); + cellIds=0; + return 0; +} + +/*! + * This method expects that 'this' is fully defined and has a spaceDim==3 and a meshDim==3. If it is not the case an exception will be thrown. + * This method returns a newly created dataarray containing cellsids in 'this' that potentially crosses the plane specified by 'origin' and 'vec'. + * @param origin is the origin of the plane. It should be an array of length 3. + * @param vec is the direction vector of the plane. It should be an array of length 3. Norm of 'vec' should be > 1e-6. + */ +DataArrayInt *MEDCouplingUMesh::getCellIdsCrossingPlane(const double *origin, const double *vec, double eps) const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=3 || getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with meshdim equal to 3 and spaceDim equal to 3 too!"); + double normm=sqrt(vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2]); + if(normm<1e-6) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellIdsCrossingPlane : parameter 'vec' should have a norm2 greater than 1e-6 !"); + double vec2[3]; + vec2[0]=vec[1]; vec2[1]=-vec[0]; vec2[2]=0.;//vec2 is the result of cross product of vec with (0,0,1) + double angle=acos(vec[2]/normm); + std::vector cellIds; + double bbox[6]; + if(angle>eps) + { + MEDCouplingAutoRefCountObjectPtr coo=_coords->deepCpy(); + MEDCouplingPointSet::Rotate3DAlg(origin,vec2,angle,coo->getNumberOfTuples(),coo->getPointer()); + MEDCouplingAutoRefCountObjectPtr mw=clone(false);//false -> shallow copy + mw->setCoords(coo); + mw->getBoundingBox(bbox); + bbox[4]=origin[2]-eps; bbox[5]=origin[2]+eps; + mw->getCellsInBoundingBox(bbox,-eps,cellIds); + } + else + { + getBoundingBox(bbox); + bbox[4]=origin[2]-eps; bbox[5]=origin[2]+eps; + getCellsInBoundingBox(bbox,-eps,cellIds); + } + MEDCouplingAutoRefCountObjectPtr ret=DataArrayInt::New(); + ret->alloc((int)cellIds.size(),1); + std::copy(cellIds.begin(),cellIds.end(),ret->getPointer()); + ret->incrRef(); + return ret; +} + /*! * This method checks that 'this' is a contiguous mesh. The user is expected to call this method on a mesh with meshdim==1. * If not an exception will thrown. If this is an empty mesh with no cell an exception will be thrown too. diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index d61087018..812ed3697 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -119,7 +119,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void renumberNodes(const int *newNodeNumbers, int newNbOfNodes); MEDCOUPLING_EXPORT void renumberNodes2(const int *newNodeNumbers, int newNbOfNodes); MEDCOUPLING_EXPORT void renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception); - MEDCOUPLING_EXPORT void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCOUPLING_EXPORT void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) const; MEDCOUPLING_EXPORT void getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox& bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureField(bool isAbs) const; MEDCOUPLING_EXPORT DataArrayDouble *getPartMeasureField(bool isAbs, const int *begin, const int *end) const; @@ -127,6 +127,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildOrthogonalField() const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildPartOrthogonalField(const int *begin, const int *end) const; MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildDirectionVectorField() const; + MEDCOUPLING_EXPORT MEDCouplingUMesh *buildSlice3D(const double *origin, const double *vec, double eps, DataArrayInt *&cellIds) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *getCellIdsCrossingPlane(const double *origin, const double *vec, double eps) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT bool isContiguous1D() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void project1D(const double *pt, const double *v, double eps, double *res) const; MEDCOUPLING_EXPORT int getCellContainingPoint(const double *pos, double eps) const; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx index 1b1b3dbb5..e0627080a 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx @@ -272,7 +272,7 @@ void MEDCouplingUMeshDesc::unserialization(const std::vector& tinyInfoD, setMeshDimension(tinyInfo[2]); } -void MEDCouplingUMeshDesc::getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) +void MEDCouplingUMeshDesc::getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) const { int dim=getSpaceDimension(); double* elem_bb=new double[2*dim]; diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx index ba62253a2..d2d8d0613 100644 --- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx +++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx @@ -65,7 +65,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT void resizeForUnserialization(const std::vector& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector& littleStrings) const; MEDCOUPLING_EXPORT void serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const; MEDCOUPLING_EXPORT void unserialization(const std::vector& tinyInfoD, const std::vector& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector& littleStrings); - MEDCOUPLING_EXPORT void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems); + MEDCOUPLING_EXPORT void getCellsInBoundingBox(const double *bbox, double eps, std::vector& elems) const; MEDCOUPLING_EXPORT void getCellsInBoundingBox(const INTERP_KERNEL::DirectedBoundingBox &bbox, double eps, std::vector& elems); MEDCOUPLING_EXPORT DataArrayInt *mergeNodes(double precision, bool& areNodesMerged, int& newNbOfNodes); MEDCOUPLING_EXPORT DataArrayInt *mergeNodes2(double precision, bool& areNodesMerged, int& newNbOfNodes); diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index 16238d612..de3651d46 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -263,6 +263,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingUMesh::findCellsIdsOnBoundary; %newobject ParaMEDMEM::MEDCouplingUMesh::getCellIdsLyingOnNodes; %newobject ParaMEDMEM::MEDCouplingUMesh::buildSetInstanceFromThis; +%newobject ParaMEDMEM::MEDCouplingUMesh::getCellIdsCrossingPlane; %newobject ParaMEDMEM::MEDCouplingUMeshCellByTypeEntry::__iter__; %newobject ParaMEDMEM::MEDCouplingUMeshCellEntry::__iter__; %newobject ParaMEDMEM::MEDCouplingExtrudedMesh::New; @@ -882,7 +883,7 @@ namespace ParaMEDMEM return ret; } - PyObject *getCellsInBoundingBox(PyObject *bbox, double eps) throw(INTERP_KERNEL::Exception) + PyObject *getCellsInBoundingBox(PyObject *bbox, double eps) const throw(INTERP_KERNEL::Exception) { std::vector elems; int size; @@ -1429,6 +1430,35 @@ namespace ParaMEDMEM PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); return ret; } + + PyObject *buildSlice3D(PyObject *origin, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception) + { + int sz; + INTERP_KERNEL::AutoPtr orig=convertPyToNewDblArr2(origin,&sz); + if(!orig || sz!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : in parameter 1 expecting origin of type list of float of size 3 !"); + INTERP_KERNEL::AutoPtr vect=convertPyToNewDblArr2(vec,&sz); + if(!vec || sz!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : in parameter 2 expecting vector of type list of float of size 3 !"); + DataArrayInt *cellIds=0; + MEDCouplingUMesh *ret0=self->buildSlice3D(orig,vect,eps,cellIds); + PyObject *ret=PyTuple_New(2); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(ret0),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(cellIds),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + return ret; + } + + DataArrayInt *getCellIdsCrossingPlane(PyObject *origin, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception) + { + int sz; + INTERP_KERNEL::AutoPtr orig=convertPyToNewDblArr2(origin,&sz); + if(!orig || sz!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellIdsCrossingPlane : in parameter 1 expecting origin of type list of float of size 3 !"); + INTERP_KERNEL::AutoPtr vect=convertPyToNewDblArr2(vec,&sz); + if(!vec || sz!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellIdsCrossingPlane : in parameter 2 expecting vector of type list of float of size 3 !"); + return self->getCellIdsCrossingPlane(orig,vect,eps); + } } void convertToPolyTypes(const std::vector& cellIdsToConvert) throw(INTERP_KERNEL::Exception); void convertAllToPoly(); -- 2.39.2