From a4dc03c6ee6d3a57d6362bdfa27767188cff0b0f Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 21 Feb 2012 14:16:30 +0000 Subject: [PATCH] buildSlice3DSurf --- src/MEDCoupling/MEDCouplingUMesh.cxx | 94 ++++++++++++- src/MEDCoupling/MEDCouplingUMesh.hxx | 1 + .../Test/MEDCouplingBasicsTest5.cxx | 132 ++++++++++++++++++ .../Test/MEDCouplingBasicsTest5.hxx | 4 + src/MEDCoupling_Swig/MEDCoupling.i | 17 +++ src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 111 +++++++++++++++ 6 files changed, 355 insertions(+), 4 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index c12ce1c5e..9707fa9aa 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2431,8 +2431,10 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const dou MEDCouplingAutoRefCountObjectPtr revDesc1=DataArrayInt::New(),revDesc2=DataArrayInt::New(); MEDCouplingAutoRefCountObjectPtr revDescIndx1=DataArrayInt::New(),revDescIndx2=DataArrayInt::New(); MEDCouplingAutoRefCountObjectPtr mDesc2=subMesh->buildDescendingConnectivity(desc2,descIndx2,revDesc2,revDescIndx2);//meshDim==2 spaceDim==3 + revDesc2=0; revDescIndx2=0; mDesc2->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds2D); MEDCouplingAutoRefCountObjectPtr mDesc1=mDesc2->buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1);//meshDim==1 spaceDim==3 + revDesc1=0; revDescIndx1=0; mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds1D); // std::vector cut3DCurve(mDesc1->getNumberOfCells(),-2); @@ -2461,7 +2463,91 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const dou } /*! - * 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 expects that 'this' is fully defined and has a spaceDim==3 and a meshDim==2. 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 different coords than 'this' and with a meshdim == 1 + * - 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 3DSurf cell id is 'this' it comes from. + * This method works only for linear meshes (non quadratic). + * If plane crosses within 'eps' a face in 'this' shared by more than 1 cell, 2 output faces will be generated. The 2 faces having the same geometry than intersecting + * face. Only 'cellIds' parameter can distinguish the 2. + * @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 3DSurf cell selection (in absolute). 'eps' is + * also used to state if new points should be created or already existing points are reused. 'eps' is also used to tells if plane overlaps a face, edge or nodes (in absolute). + */ +MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3DSurf(const double *origin, const double *vec, double eps, DataArrayInt *&cellIds) const throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=2 || getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with meshdim equal to 2 and spaceDim equal to 3 !"); + MEDCouplingAutoRefCountObjectPtr candidates=getCellIdsCrossingPlane(origin,vec,eps); + if(candidates->empty()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D : No 3D cells in this intercepts the specified plane considering bounding boxes !"); + std::vector nodes; + std::vector cellIds1D; + MEDCouplingAutoRefCountObjectPtr subMesh=static_cast(buildPartOfMySelf(candidates->begin(),candidates->end(),false)); + int nbNodes1=subMesh->getNumberOfNodes(); + subMesh->findNodesOnPlane(origin,vec,eps,nodes); + MEDCouplingAutoRefCountObjectPtr desc1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr descIndx1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr revDesc1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr revDescIndx1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr mDesc1=subMesh->buildDescendingConnectivity(desc1,descIndx1,revDesc1,revDescIndx1);//meshDim==1 spaceDim==3 + mDesc1->fillCellIdsToKeepFromNodeIds(&nodes[0],&nodes[0]+nodes.size(),true,cellIds1D); + // + std::vector cut3DCurve(mDesc1->getNumberOfCells(),-2); + for(std::vector::const_iterator it=cellIds1D.begin();it!=cellIds1D.end();it++) + cut3DCurve[*it]=-1; + mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve); + int ncellsSub=subMesh->getNumberOfCells(); + std::vector< std::pair > cut3DSurf(ncellsSub); + AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,subMesh->getNodalConnectivity()->getConstPointer(),subMesh->getNodalConnectivityIndex()->getConstPointer(), + mDesc1->getNodalConnectivity()->getConstPointer(),mDesc1->getNodalConnectivityIndex()->getConstPointer(), + desc1->getConstPointer(),descIndx1->getConstPointer(),cut3DSurf); + std::vector conn,connI,cellIds2; connI.push_back(0); + const int *nodal=subMesh->getNodalConnectivity()->getConstPointer(); + const int *nodalI=subMesh->getNodalConnectivityIndex()->getConstPointer(); + for(int i=0;i ret=MEDCouplingUMesh::New("Slice3DSurf",1); + ret->setCoords(mDesc1->getCoords()); + MEDCouplingAutoRefCountObjectPtr c=DataArrayInt::New(); + c->alloc((int)conn.size(),1); std::copy(conn.begin(),conn.end(),c->getPointer()); + MEDCouplingAutoRefCountObjectPtr cI=DataArrayInt::New(); + cI->alloc((int)connI.size(),1); std::copy(connI.begin(),connI.end(),cI->getPointer()); + ret->setConnectivity(c,cI,true); + cellIds=candidates->selectByTupleId(&cellIds2[0],&cellIds2[0]+cellIds2.size()); + ret->incrRef(); + return ret; +} + +/*! + * This method expects that 'this' is fully defined and has a spaceDim==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. @@ -2469,8 +2555,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildSlice3D(const double *origin, const dou 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!"); + if(getSpaceDimension()!=3) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildSlice3D works on umeshes with spaceDim equal to 3 !"); 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 !"); @@ -5447,7 +5533,7 @@ void MEDCouplingUMesh::AssemblyForSplitFrom3DCurve(const std::vector& cut3D {// case when plane is on a multi colinear edge of a polyhedron if(res.size()==2*nbOfSeg) { - cut3DSurf[i].first=-2; cut3DSurf[i].first=i; + cut3DSurf[i].first=-2; cut3DSurf[i].second=i; } else throw INTERP_KERNEL::Exception("MEDCouplingUMesh::AssemblyPointsFrom3DCurve : unexpected situation !"); diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 0c8a1cd40..4a124842e 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -128,6 +128,7 @@ namespace ParaMEDMEM 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 MEDCouplingUMesh *buildSlice3DSurf(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; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx index 3a8656e3d..c49d78d8f 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.cxx @@ -172,3 +172,135 @@ void MEDCouplingBasicsTest5::testGetCellIdsCrossingPlane1() mesh2D->decrRef(); } +void MEDCouplingBasicsTest5::testBuildSlice3D1() +{ + MEDCouplingUMesh *mesh2D=0; + MEDCouplingUMesh *mesh3D=build3DExtrudedUMesh_1(mesh2D); + mesh2D->decrRef(); + // First slice in the middle of 3D cells + const double vec1[3]={-0.07,1.,0.07}; + const double origin1[3]={1.524,1.4552,1.74768}; + DataArrayInt *ids=0; + MEDCouplingUMesh *slice1=mesh3D->buildSlice3D(origin1,vec1,1e-10,ids); + const int expected1[9]={1,3,4,7,9,10,13,15,16}; + const int expected2[47]={5,42,41,40,43,44,5,42,46,45,41,5,44,43,40,47,48,5,49,42,44,50,5,49,51,46,42,5,50,44,48,52,5,53,49,50,54,5,53,55,51,49,5,54,50,52,56}; + const int expected3[10]={0,6,11,17,22,27,32,37,42,47}; + const double expected4[171]={1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,2.,2.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,2.,2.,1.,1.,1.,2.,1.,1.25,2.,1.,1.5,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,1.,1.,3.,1.,1.25,3.,1.,1.5,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,1.,1.5408576,0.,2.,1.6108576000000001,0.,2.,1.5408576,1.,1.,1.5,0.5836800000000008,1.,1.4708576,1.,3.,1.6808576,0.,3.,1.6108576000000001,1.,0.,1.4708576,0.,0.,1.4008576,1.,2.,1.4708576,2.,1.,1.4008576000000001,2.,3.,1.5408575999999998,2.,0.,1.3308575999999999,2.,2.,1.4008576,3.,1.,1.3308576,3.,3.,1.4708576,3.,0.,1.2608576,3.}; + CPPUNIT_ASSERT_EQUAL(2,slice1->getMeshDimension()); + CPPUNIT_ASSERT_EQUAL(3,slice1->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(57,slice1->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(9,slice1->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(9,ids->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(47,slice1->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(10,slice1->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expected1,expected1+9,ids->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected2,expected2+47,slice1->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected3,expected3+10,slice1->getNodalConnectivityIndex()->getConstPointer())); + for(int i=0;i<171;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected4[i],slice1->getCoords()->getIJ(0,i),1e-12); + ids->decrRef(); + slice1->decrRef(); + // 2nd slice based on already existing nodes of mesh3D. + const double vec2[3]={0.,3.,1.}; + const double origin2[3]={2.5,1.,3.}; + slice1=mesh3D->buildSlice3D(origin2,vec2,1e-10,ids); + const int expected5[49]={5,50,10,4,51,5,50,52,7,10,5,51,4,5,53,5,54,50,51,55,56,5,54,57,52,50,5,56,55,51,53,58,5,38,59,56,54,43,5,54,57,46,43,5,38,59,56,58,48}; + const int expected6[10]={0,5,10,15,21,26,32,38,43,49}; + const double expected7[180]={1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,1.,3.,0.,2.,2.,0.,2.,3.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25,2.,1.,0.,2.,1.,1.5,2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,0.,0.,3.,1.,1.,3.,1.,1.25,3.,1.,0.,3.,1.,1.5,3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,2.,1.6666666666666667,1.,1.,1.6666666666666667,1.,3.,1.6666666666666667,1.,0.,1.6666666666666667,1.,2.,1.3333333333333335,2.,1.,1.5,1.5,1.,1.3333333333333333,2.,3.,1.3333333333333335,2.,0.,1.3333333333333335,2.,1.,1.25,2.25}; + CPPUNIT_ASSERT_EQUAL(2,slice1->getMeshDimension()); + CPPUNIT_ASSERT_EQUAL(3,slice1->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(60,slice1->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(9,slice1->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(9,ids->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(49,slice1->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(10,slice1->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expected1,expected1+9,ids->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected5,expected5+49,slice1->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected6,expected6+10,slice1->getNodalConnectivityIndex()->getConstPointer())); + for(int i=0;i<180;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected7[i],slice1->getCoords()->getIJ(0,i),1e-12); + ids->decrRef(); + slice1->decrRef(); + // 3rd slice based on shared face of mesh3D. + const double vec3[3]={0.,0.,1.}; + const double origin3[3]={2.5,1.,2.}; + slice1=mesh3D->buildSlice3D(origin3,vec3,1e-10,ids); + const int expected8[12]={6,7,8,9,10,11,12,13,14,15,16,17}; + const int expected9[68]={5,15,26,16,18,5,16,21,28,22,19,17,5,18,20,21,16,5,21,24,25,28,5,26,16,17,19,22,23,5,22,27,29,28,5,15,26,16,18,5,16,21,28,22,19,17,5,18,20,21,16,5,21,24,25,28,5,26,16,17,19,22,23,5,22,27,29,28}; + const int expected10[13]={0,5,12,17,22,29,34,39,46,51,56,63,68}; + const double expected11[135]={0.,0.,1.,1.,1.,1.,1.,1.25, 1.,1.,0.,1.,1.,1.5, 1.,2.,0.,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25, 2.,1.,0.,2.,1.,1.5, 2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,1.,3.,2.,2.,2.,2.,2.,3.,2.,0.,0.,3.,1.,1.,3.,1.,1.25, 3.,1.,0.,3.,1.,1.5, 3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,1.,3.,3.,2.,2.,3.,2.,3.,3.}; + CPPUNIT_ASSERT_EQUAL(2,slice1->getMeshDimension()); + CPPUNIT_ASSERT_EQUAL(3,slice1->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(45,slice1->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(12,slice1->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(12,ids->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(68,slice1->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(13,slice1->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expected8,expected8+12,ids->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected9,expected9+68,slice1->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected10,expected10+13,slice1->getNodalConnectivityIndex()->getConstPointer())); + for(int i=0;i<135;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected11[i],slice1->getCoords()->getIJ(0,i),1e-12); + ids->decrRef(); + slice1->decrRef(); + // + mesh3D->decrRef(); +} + +void MEDCouplingBasicsTest5::testBuildSlice3DSurf1() +{ + MEDCouplingUMesh *mesh2D=0; + MEDCouplingUMesh *mesh3D=build3DExtrudedUMesh_1(mesh2D); + mesh2D->decrRef(); + DataArrayInt *a=DataArrayInt::New(),*b=DataArrayInt::New(),*c=DataArrayInt::New(),*d=DataArrayInt::New(); + mesh2D=mesh3D->buildDescendingConnectivity(a,b,c,d); + a->decrRef(); b->decrRef(); c->decrRef(); d->decrRef(); + mesh3D->decrRef(); + // + const double vec1[3]={-0.07,1.,0.07}; + const double origin1[3]={1.524,1.4552,1.74768}; + DataArrayInt *ids=0; + MEDCouplingUMesh *slice1=mesh2D->buildSlice3DSurf(origin1,vec1,1e-10,ids); + const int expected1[25]={6,8,10,11,13,18,19,21,23,25,26,38,41,43,47,49,52,53,64,67,69,73,75,78,79}; + const int expected2[75]={1,40,41,1,42,41,1,40,43,1,44,43,1,42,44,1,45,41,1,42,46,1,46,45,1,47,40,1,47,48,1,44,48,1,49,42,1,44,50,1,49,50,1,49,51,1,51,46,1,48,52,1,50,52,1,53,49,1,50,54,1,53,54,1,53,55,1,55,51,1,52,56,1,54,56}; + const int expected3[26]={0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60,63,66,69,72,75}; + const double expected4[171]={1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,2.,2.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,2.,2.,1.,1.,1.,2.,1.,1.25,2.,1.,1.5,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,1.,1.,3.,1.,1.25,3.,1.,1.5,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,1.,1.5408576,0.,2.,1.6108576000000001,0.,2.,1.5408576,1.,1.,1.5,0.5836800000000008,1.,1.4708576,1.,3.,1.6808576,0.,3.,1.6108576000000001,1.,0.,1.4708576,0.,0.,1.4008576,1.,2.,1.4708576,2.,1.,1.4008576000000001,2.,3.,1.5408575999999998,2.,0.,1.3308575999999999,2.,2.,1.4008576,3.,1.,1.3308576,3.,3.,1.4708576,3.,0.,1.2608576,3.}; + CPPUNIT_ASSERT_EQUAL(1,slice1->getMeshDimension()); + CPPUNIT_ASSERT_EQUAL(3,slice1->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(57,slice1->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(25,slice1->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(25,ids->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(75,slice1->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(26,slice1->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expected1,expected1+25,ids->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected2,expected2+47,slice1->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected3,expected3+26,slice1->getNodalConnectivityIndex()->getConstPointer())); + for(int i=0;i<171;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected4[i],slice1->getCoords()->getIJ(0,i),1e-12); + ids->decrRef(); + slice1->decrRef(); + // + const double vec2[3]={0.,0.,1.}; + const double origin2[3]={2.5,1.,2.}; + slice1=mesh2D->buildSlice3DSurf(origin2,vec2,1e-10,ids); + const int expected5[68]={32,32,32,32,33,34,35,36,37,38,39,40,41,42,43,43,43,43,43,43,44,44,44,44,45,46,47,47,47,47,48,49,50,51,52,53,53,53,53,53,53,54,54,54,54,55,56,57,59,60,61,62,63,64,65,66,67,68,71,72,74,75,76,77,78,81,82,83}; + const int expected6[204]={1,15,18,1,18,16,1,16,26,1,26,15,1,26,15,1,16,26,1,18,16,1,15,18,1,16,21,1,21,28,1,22,28,1,19,22,1,17,19,1,16,17,1,16,21,1,21,28,1,28,22,1,22,19,1,19,17,1,17,16,1,16,18,1,18,20,1,20,21,1,21,16,1,20,21,1,18,20,1,28,21,1,21,24,1,24,25,1,25,28,1,25,28,1,24,25,1,21,24,1,23,22,1,26,23,1,26,16,1,16,17,1,17,19,1,19,22,1,22,23,1,23,26,1,22,28,1,28,29,1,29,27,1,27,22,1,27,22,1,29,27,1,28,29,1,26,15,1,16,26,1,18,16,1,15,18,1,16,21,1,21,28,1,22,28,1,19,22,1,17,19,1,16,17,1,20,21,1,18,20,1,25,28,1,24,25,1,21,24,1,23,22,1,26,23,1,27,22,1,29,27,1,28,29}; + const int expected7[69]={0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60,63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114,117,120,123,126,129,132,135,138,141,144,147,150,153,156,159,162,165,168,171,174,177,180,183,186,189,192,195,198,201,204}; + const double expected8[135]={0.,0.,1.,1.,1.,1.,1.,1.25, 1.,1.,0.,1.,1.,1.5, 1.,2.,0.,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25, 2.,1.,0.,2.,1.,1.5, 2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,1.,3.,2.,2.,2.,2.,2.,3.,2.,0.,0.,3.,1.,1.,3.,1.,1.25, 3.,1.,0.,3.,1.,1.5, 3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,1.,3.,3.,2.,2.,3.,2.,3.,3.}; + CPPUNIT_ASSERT_EQUAL(1,slice1->getMeshDimension()); + CPPUNIT_ASSERT_EQUAL(3,slice1->getSpaceDimension()); + CPPUNIT_ASSERT_EQUAL(45,slice1->getNumberOfNodes()); + CPPUNIT_ASSERT_EQUAL(68,slice1->getNumberOfCells()); + CPPUNIT_ASSERT_EQUAL(68,ids->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(204,slice1->getNodalConnectivity()->getNumberOfTuples()); + CPPUNIT_ASSERT_EQUAL(69,slice1->getNodalConnectivityIndex()->getNumberOfTuples()); + CPPUNIT_ASSERT(std::equal(expected5,expected5+68,ids->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected6,expected6+171,slice1->getNodalConnectivity()->getConstPointer())); + CPPUNIT_ASSERT(std::equal(expected7,expected7+69,slice1->getNodalConnectivityIndex()->getConstPointer())); + for(int i=0;i<135;i++) + CPPUNIT_ASSERT_DOUBLES_EQUAL(expected8[i],slice1->getCoords()->getIJ(0,i),1e-12); + ids->decrRef(); + slice1->decrRef(); + // + mesh2D->decrRef(); +} diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx index 6ee46d0c6..6d8a20610 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest5.hxx @@ -38,11 +38,15 @@ namespace ParaMEDMEM CPPUNIT_TEST( testUMeshTessellate2D1 ); CPPUNIT_TEST( testIntersect2DMeshesTmp4 ); CPPUNIT_TEST( testGetCellIdsCrossingPlane1 ); + CPPUNIT_TEST( testBuildSlice3D1 ); + CPPUNIT_TEST( testBuildSlice3DSurf1 ); CPPUNIT_TEST_SUITE_END(); public: void testUMeshTessellate2D1(); void testIntersect2DMeshesTmp4(); void testGetCellIdsCrossingPlane1(); + void testBuildSlice3D1(); + void testBuildSlice3DSurf1(); }; } diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index de3651d46..fc75b9d39 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -1448,6 +1448,23 @@ namespace ParaMEDMEM return ret; } + PyObject *buildSlice3DSurf(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->buildSlice3DSurf(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; diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index a1988d687..bb038c49b 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -8604,6 +8604,117 @@ class MEDCouplingBasicsTest(unittest.TestCase): ids2=mesh3D.getCellIdsCrossingPlane(origin,vec2,1e-10) self.assertEqual([6,7,8,9,10,11],ids2.getValues()) pass + + def testBuildSlice3D1(self): + mesh3D,mesh2D=MEDCouplingDataForTest.build3DExtrudedUMesh_1(); + vec1=[-0.07,1.,0.07] + origin1=[1.524,1.4552,1.74768] + slice1,ids=mesh3D.buildSlice3D(origin1,vec1,1e-10); + expected1=[1,3,4,7,9,10,13,15,16] + expected2=[5,42,41,40,43,44,5,42,46,45,41,5,44,43,40,47,48,5,49,42,44,50,5,49,51,46,42,5,50,44,48,52,5,53,49,50,54,5,53,55,51,49,5,54,50,52,56] + expected3=[0,6,11,17,22,27,32,37,42,47] + expected4=[1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,2.,2.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,2.,2.,1.,1.,1.,2.,1.,1.25,2.,1.,1.5,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,1.,1.,3.,1.,1.25,3.,1.,1.5,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,1.,1.5408576,0.,2.,1.6108576000000001,0.,2.,1.5408576,1.,1.,1.5,0.5836800000000008,1.,1.4708576,1.,3.,1.6808576,0.,3.,1.6108576000000001,1.,0.,1.4708576,0.,0.,1.4008576,1.,2.,1.4708576,2.,1.,1.4008576000000001,2.,3.,1.5408575999999998,2.,0.,1.3308575999999999,2.,2.,1.4008576,3.,1.,1.3308576,3.,3.,1.4708576,3.,0.,1.2608576,3.] + self.assertEqual(2,slice1.getMeshDimension()); + self.assertEqual(3,slice1.getSpaceDimension()); + self.assertEqual(57,slice1.getNumberOfNodes()); + self.assertEqual(9,slice1.getNumberOfCells()); + self.assertEqual(9,ids.getNumberOfTuples()); + self.assertEqual(47,slice1.getNodalConnectivity().getNumberOfTuples()); + self.assertEqual(10,slice1.getNodalConnectivityIndex().getNumberOfTuples()); + self.assertEqual(expected1,ids.getValues()); + self.assertEqual(expected2,slice1.getNodalConnectivity().getValues()); + self.assertEqual(expected3,slice1.getNodalConnectivityIndex().getValues()); + for i in xrange(171): + self.assertAlmostEqual(expected4[i],slice1.getCoords().getIJ(0,i),12); + pass + # 2nd slice based on already existing nodes of mesh3D. + vec2=[0.,3.,1.] + origin2=[2.5,1.,3.] + slice1,ids=mesh3D.buildSlice3D(origin2,vec2,1e-10); + expected5=[5,50,10,4,51,5,50,52,7,10,5,51,4,5,53,5,54,50,51,55,56,5,54,57,52,50,5,56,55,51,53,58,5,38,59,56,54,43,5,54,57,46,43,5,38,59,56,58,48] + expected6=[0,5,10,15,21,26,32,38,43,49] + expected7=[1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,1.,3.,0.,2.,2.,0.,2.,3.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25,2.,1.,0.,2.,1.,1.5,2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,0.,0.,3.,1.,1.,3.,1.,1.25,3.,1.,0.,3.,1.,1.5,3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,2.,1.6666666666666667,1.,1.,1.6666666666666667,1.,3.,1.6666666666666667,1.,0.,1.6666666666666667,1.,2.,1.3333333333333335,2.,1.,1.5,1.5,1.,1.3333333333333333,2.,3.,1.3333333333333335,2.,0.,1.3333333333333335,2.,1.,1.25,2.25] + self.assertEqual(2,slice1.getMeshDimension()); + self.assertEqual(3,slice1.getSpaceDimension()); + self.assertEqual(60,slice1.getNumberOfNodes()); + self.assertEqual(9,slice1.getNumberOfCells()); + self.assertEqual(9,ids.getNumberOfTuples()); + self.assertEqual(49,slice1.getNodalConnectivity().getNumberOfTuples()); + self.assertEqual(10,slice1.getNodalConnectivityIndex().getNumberOfTuples()); + self.assertEqual(expected1,ids.getValues()); + self.assertEqual(expected5,slice1.getNodalConnectivity().getValues()); + self.assertEqual(expected6,slice1.getNodalConnectivityIndex().getValues()); + for i in xrange(180): + self.assertAlmostEqual(expected7[i],slice1.getCoords().getIJ(0,i),12); + pass + # 3rd slice based on shared face of mesh3D. + vec3=[0.,0.,1.] + origin3=[2.5,1.,2.] + slice1,ids=mesh3D.buildSlice3D(origin3,vec3,1e-10); + expected8=[6,7,8,9,10,11,12,13,14,15,16,17] + expected9=[5,15,26,16,18,5,16,21,28,22,19,17,5,18,20,21,16,5,21,24,25,28,5,26,16,17,19,22,23,5,22,27,29,28,5,15,26,16,18,5,16,21,28,22,19,17,5,18,20,21,16,5,21,24,25,28,5,26,16,17,19,22,23,5,22,27,29,28] + expected10=[0,5,12,17,22,29,34,39,46,51,56,63,68] + expected11=[0.,0.,1.,1.,1.,1.,1.,1.25,1.,1.,0.,1.,1.,1.5,1.,2.,0.,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25,2.,1.,0.,2.,1.,1.5,2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,1.,3.,2.,2.,2.,2.,2.,3.,2.,0.,0.,3.,1.,1.,3.,1.,1.25,3.,1.,0.,3.,1.,1.5,3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,1.,3.,3.,2.,2.,3.,2.,3.,3.] + self.assertEqual(2,slice1.getMeshDimension()); + self.assertEqual(3,slice1.getSpaceDimension()); + self.assertEqual(45,slice1.getNumberOfNodes()); + self.assertEqual(12,slice1.getNumberOfCells()); + self.assertEqual(12,ids.getNumberOfTuples()); + self.assertEqual(68,slice1.getNodalConnectivity().getNumberOfTuples()); + self.assertEqual(13,slice1.getNodalConnectivityIndex().getNumberOfTuples()); + self.assertEqual(expected8,ids.getValues()); + self.assertEqual(expected9,slice1.getNodalConnectivity().getValues()); + self.assertEqual(expected10,slice1.getNodalConnectivityIndex().getValues()); + for i in xrange(135): + self.assertAlmostEqual(expected11[i],slice1.getCoords().getIJ(0,i),12); + pass + pass + + def testBuildSlice3DSurf1(self): + mesh3D,mesh2D=MEDCouplingDataForTest.build3DExtrudedUMesh_1(); + mesh2D=mesh3D.buildDescendingConnectivity()[0]; + vec1=[-0.07,1.,0.07] + origin1=[1.524,1.4552,1.74768] + slice1,ids=mesh2D.buildSlice3DSurf(origin1,vec1,1e-10); + expected1=[6,8,10,11,13,18,19,21,23,25,26,38,41,43,47,49,52,53,64,67,69,73,75,78,79] + expected2=[1,40,41,1,42,41,1,40,43,1,44,43,1,42,44,1,45,41,1,42,46,1,46,45,1,47,40,1,47,48,1,44,48,1,49,42,1,44,50,1,49,50,1,49,51,1,51,46,1,48,52,1,50,52,1,53,49,1,50,54,1,53,54,1,53,55,1,55,51,1,52,56,1,54,56] + expected3=[0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60,63,66,69,72,75]; + expected4=[1.,1.,0.,1.,1.25,0.,1.,1.5,0.,2.,1.,0.,1.,2.,0.,0.,2.,0.,3.,1.,0.,3.,2.,0.,0.,1.,0.,2.,2.,0.,1.,1.,1.,1.,1.25,1.,1.,1.5,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,2.,2.,1.,1.,1.,2.,1.,1.25,2.,1.,1.5,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,2.,2.,2.,1.,1.,3.,1.,1.25,3.,1.,1.5,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,2.,2.,3.,1.,1.5408576,0.,2.,1.6108576000000001,0.,2.,1.5408576,1.,1.,1.5,0.5836800000000008,1.,1.4708576,1.,3.,1.6808576,0.,3.,1.6108576000000001,1.,0.,1.4708576,0.,0.,1.4008576,1.,2.,1.4708576,2.,1.,1.4008576000000001,2.,3.,1.5408575999999998,2.,0.,1.3308575999999999,2.,2.,1.4008576,3.,1.,1.3308576,3.,3.,1.4708576,3.,0.,1.2608576,3.] + self.assertEqual(1,slice1.getMeshDimension()); + self.assertEqual(3,slice1.getSpaceDimension()); + self.assertEqual(57,slice1.getNumberOfNodes()); + self.assertEqual(25,slice1.getNumberOfCells()); + self.assertEqual(25,ids.getNumberOfTuples()); + self.assertEqual(75,slice1.getNodalConnectivity().getNumberOfTuples()); + self.assertEqual(26,slice1.getNodalConnectivityIndex().getNumberOfTuples()); + self.assertEqual(expected1,ids.getValues()); + self.assertEqual(expected2,slice1.getNodalConnectivity().getValues()); + self.assertEqual(expected3,slice1.getNodalConnectivityIndex().getValues()); + for i in xrange(171): + self.assertAlmostEqual(expected4[i],slice1.getCoords().getIJ(0,i),12); + pass + # + vec2=[0.,0.,1.] + origin2=[2.5,1.,2.] + slice1,ids=mesh2D.buildSlice3DSurf(origin2,vec2,1e-10); + expected5=[32,32,32,32,33,34,35,36,37,38,39,40,41,42,43,43,43,43,43,43,44,44,44,44,45,46,47,47,47,47,48,49,50,51,52,53,53,53,53,53,53,54,54,54,54,55,56,57,59,60,61,62,63,64,65,66,67,68,71,72,74,75,76,77,78,81,82,83] + expected6=[1,15,18,1,18,16,1,16,26,1,26,15,1,26,15,1,16,26,1,18,16,1,15,18,1,16,21,1,21,28,1,22,28,1,19,22,1,17,19,1,16,17,1,16,21,1,21,28,1,28,22,1,22,19,1,19,17,1,17,16,1,16,18,1,18,20,1,20,21,1,21,16,1,20,21,1,18,20,1,28,21,1,21,24,1,24,25,1,25,28,1,25,28,1,24,25,1,21,24,1,23,22,1,26,23,1,26,16,1,16,17,1,17,19,1,19,22,1,22,23,1,23,26,1,22,28,1,28,29,1,29,27,1,27,22,1,27,22,1,29,27,1,28,29,1,26,15,1,16,26,1,18,16,1,15,18,1,16,21,1,21,28,1,22,28,1,19,22,1,17,19,1,16,17,1,20,21,1,18,20,1,25,28,1,24,25,1,21,24,1,23,22,1,26,23,1,27,22,1,29,27,1,28,29] + expected7=[0,3,6,9,12,15,18,21,24,27,30,33,36,39,42,45,48,51,54,57,60,63,66,69,72,75,78,81,84,87,90,93,96,99,102,105,108,111,114,117,120,123,126,129,132,135,138,141,144,147,150,153,156,159,162,165,168,171,174,177,180,183,186,189,192,195,198,201,204]; + expected8=[0.,0.,1.,1.,1.,1.,1.,1.25, 1.,1.,0.,1.,1.,1.5, 1.,2.,0.,1.,2.,1.,1.,1.,2.,1.,0.,2.,1.,3.,1.,1.,3.,2.,1.,0.,1.,1.,1.,3.,1.,2.,2.,1.,2.,3.,1.,0.,0.,2.,1.,1.,2.,1.,1.25, 2.,1.,0.,2.,1.,1.5, 2.,2.,0.,2.,2.,1.,2.,1.,2.,2.,0.,2.,2.,3.,1.,2.,3.,2.,2.,0.,1.,2.,1.,3.,2.,2.,2.,2.,2.,3.,2.,0.,0.,3.,1.,1.,3.,1.,1.25, 3.,1.,0.,3.,1.,1.5, 3.,2.,0.,3.,2.,1.,3.,1.,2.,3.,0.,2.,3.,3.,1.,3.,3.,2.,3.,0.,1.,3.,1.,3.,3.,2.,2.,3.,2.,3.,3.] + self.assertEqual(1,slice1.getMeshDimension()); + self.assertEqual(3,slice1.getSpaceDimension()); + self.assertEqual(45,slice1.getNumberOfNodes()); + self.assertEqual(68,slice1.getNumberOfCells()); + self.assertEqual(68,ids.getNumberOfTuples()); + self.assertEqual(204,slice1.getNodalConnectivity().getNumberOfTuples()); + self.assertEqual(69,slice1.getNodalConnectivityIndex().getNumberOfTuples()); + self.assertEqual(expected5,ids.getValues()); + self.assertEqual(expected6,slice1.getNodalConnectivity().getValues()); + self.assertEqual(expected7,slice1.getNodalConnectivityIndex().getValues()); + for i in xrange(135): + self.assertAlmostEqual(expected8[i],slice1.getCoords().getIJ(0,i),12); + pass + pass def setUp(self): pass -- 2.39.2