From 745be713d23f563c12a6e695b1686a22d53f833c Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 13 Feb 2012 16:10:19 +0000 Subject: [PATCH] IntersectMesh with quadratic elements --- .../Geometric2D/InterpKernelGeo2DEdge.hxx | 4 +- .../InterpKernelGeo2DEdgeArcCircle.cxx | 6 +- .../InterpKernelGeo2DElementaryEdge.cxx | 2 +- .../InterpKernelGeo2DElementaryEdge.hxx | 2 +- .../InterpKernelGeo2DQuadraticPolygon.cxx | 80 +++++++++++++++---- .../InterpKernelGeo2DQuadraticPolygon.hxx | 3 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 22 +++-- src/MEDCoupling/MEDCouplingUMesh.hxx | 2 +- src/MEDLoader/MEDFileMesh.cxx | 8 ++ src/MEDLoader/MEDFileMesh.hxx | 1 + src/MEDLoader/Swig/MEDLoader.i | 1 + 11 files changed, 103 insertions(+), 28 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index e712a0d0d..cdd81b2c2 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -261,7 +261,8 @@ namespace INTERP_KERNEL std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const = 0; virtual void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const = 0; - virtual void sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis) = 0; + virtual void sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis) = 0; + virtual Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const = 0; protected: Edge():_cnt(1),_loc(FULL_UNKNOWN),_start(0),_end(0) { } virtual ~Edge(); @@ -271,7 +272,6 @@ namespace INTERP_KERNEL //! The code 'code' is built by method combineCodes static bool SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *nE, bool direction, int code, ComposedEdge& outVal1, ComposedEdge& outVal2); - virtual Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const = 0; protected: mutable unsigned char _cnt; mutable TypeOfEdgeLocInPolygon _loc; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index b9a8f0cdd..640f5b122 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -497,18 +497,18 @@ void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int res stream << "5 1 0 1 "; fillXfigStreamForLoc(stream); stream << " 7 50 -1 -1 0.000 0 "; - if( (direction && _angle>=0) || (!direction && _angle<0)) + if( (direction && (-_angle)>=0) || (!direction && (-_angle)<0)) stream << '0';//'0' else stream << '1';//'1' - stream << " 0 0 "; + stream << " 1 0 "; stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " "; direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box); Node *middle=buildRepresentantOfMySelf(); middle->dumpInXfigFile(stream,resolution,box); middle->decrRef(); direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box); - stream << std::endl; + stream << std::endl << "1 1 2.00 120.00 180.00" << std::endl; } void EdgeArcCircle::update(Node *m) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx index 0294997f2..59100475a 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx @@ -224,7 +224,7 @@ void ElementaryEdge::fillGlobalInfoAbs2(const std::map& edgesThis, std::vector& addCoo, std::map mapAddCoo) const; void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, std::vector& edgesOther, std::vector& addCoo, std::map& mapAddCoo) const; - static ElementaryEdge *BuildEdgeFromCrudeDataArray(bool isQuad, bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end); + static ElementaryEdge *BuildEdgeFromCrudeDataArray(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end); private: bool _direction; Edge *_ptr; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 178a184ac..1a96efa8e 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -292,21 +292,73 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, const std::map& mapp, bool isQuad, const int *descBg, const int *descEnd, const std::vector >& intersectEdges) +void QuadraticPolygon::buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, + const int *descBg, const int *descEnd, const std::vector >& intersectEdges) { - std::size_t nbOfSeg=std::distance(descBg,descEnd); - for(std::size_t i=0;i0; - int edgeId=abs(descBg[i])-1; - const std::vector& subEdge=intersectEdges[edgeId]; - std::size_t nbOfSubEdges=subEdge.size()/2; - for(std::size_t j=0;j0; + int edgeId=abs(descBg[i])-1; + const std::vector& subEdge=intersectEdges[edgeId]; + std::size_t nbOfSubEdges=subEdge.size()/2; + for(std::size_t j=0;j0; + int edgeId=abs(descBg[i])-1; + const std::vector& subEdge=intersectEdges[edgeId]; + std::size_t nbOfSubEdges=subEdge.size()/2; + if(colinearity) + { + for(std::size_t j=0;jbuildEdgeLyingOnMe(start,end); + ElementaryEdge *eee=new ElementaryEdge(ee,true); + pushBack(eee); + } + e->decrRef(); + } + st0->decrRef(); endd0->decrRef(); middle0->decrRef(); } } } @@ -350,7 +402,7 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map build new Edge instance Node *start=(*mapp.find(idBg)).second; Node *end=(*mapp.find(idEnd)).second; - ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(isQuad,true,start,end); + ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end); pushBack(e); } else diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index 8ccbde27f..bab62ab32 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -56,7 +56,8 @@ namespace INTERP_KERNEL double intersectWithAbs(QuadraticPolygon& other, double* barycenter); void splitAbs(QuadraticPolygon& other, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, const std::vector& otherEdgeIds, std::vector& edgesThis, int cellIdThis, std::vector< std::vector >& edgesInOtherColinearWithThis, std::vector< std::vector >& subDivOther, std::vector& addCoo); - void buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const int *descBg, const int *descEnd, const std::vector >& intersectEdges); + void buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const int *nodalBg, const double *coords, + const int *descBg, const int *descEnd, const std::vector >& intersectEdges); void buildFromCrudeDataArray2(const std::map& mapp, bool isQuad, const int *descBg, const int *descEnd, const std::vector >& intersectEdges, const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector >& intersectEdges1, const std::vector< std::vector >& colinear1); diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 7e920f6ad..9a517ec51 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2527,7 +2527,15 @@ namespace ParaMEDMEM } case INTERP_KERNEL::NORM_SEG3: { - ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]],mapp2[bg[2]],mapp2[bg[1]]); + INTERP_KERNEL::EdgeLin * e1=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]],mapp2[bg[2]]); + INTERP_KERNEL::EdgeLin *e2=new INTERP_KERNEL::EdgeLin(mapp2[bg[2]],mapp2[bg[1]]); + INTERP_KERNEL::SegSegIntersector inters(*e1,*e2); + bool colinearity=inters.areColinears(); + delete e1; delete e2; + if(colinearity) + ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]],mapp2[bg[1]]); + else + ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]],mapp2[bg[2]],mapp2[bg[1]]); break; } default: @@ -4825,12 +4833,11 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 MEDCouplingAutoRefCountObjectPtr dd5(m1Desc),dd6(m2Desc); std::vector< std::vector > intersectEdge2; BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2); - std::vector b1=m1Desc->getQuadraticStatus(); std::vector b2=m1Desc->getQuadraticStatus(); subDiv2.clear(); dd5=0; dd6=0; std::vector cr,crI; std::vector cNb1,cNb2; - BuildIntersecting2DCellsFromEdges(eps,m1,b1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,colinear2,m2,b2,desc2->getConstPointer(),descIndx2->getConstPointer(),intersectEdge2,addCoo, + BuildIntersecting2DCellsFromEdges(eps,m1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,colinear2,m2,b2,desc2->getConstPointer(),descIndx2->getConstPointer(),intersectEdge2,addCoo, /* outputs -> */cr,crI,cNb1,cNb2); // MEDCouplingAutoRefCountObjectPtr addCooDa=DataArrayDouble::New(); @@ -4852,7 +4859,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1 /// @endcond -void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const std::vector& b1, const int *desc1, const int *descIndx1, +void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, const MEDCouplingUMesh *m2, const std::vector& b2, const int *desc2, const int *descIndx2, const std::vector >& intesctEdges2, const std::vector& addCoords, @@ -4861,6 +4868,8 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo static const int SPACEDIM=2; std::vector bbox1,bbox2; const double *coo1=m1->getCoords()->getConstPointer(); + const int *conn1=m1->getNodalConnectivity()->getConstPointer(); + const int *connI1=m1->getNodalConnectivityIndex()->getConstPointer(); int offset1=m1->getNumberOfNodes(); const double *coo2=m2->getCoords()->getConstPointer(); int offset2=offset1+m2->getNumberOfNodes(); @@ -4876,8 +4885,11 @@ void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCo std::map mapp; std::map mappRev; INTERP_KERNEL::QuadraticPolygon pol1; + INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]]; + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ); MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev); - pol1.buildFromCrudeDataArray(mappRev,b1[i],desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); + pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1, + desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1); std::vector crTmp,crITmp; crITmp.push_back(crI.back()); for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++) diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index d9aec073c..dc84121d8 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -220,7 +220,7 @@ namespace ParaMEDMEM MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2, std::vector& addCoo) throw(INTERP_KERNEL::Exception); static void BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, const std::vector& addCoo, const std::vector< std::vector >& subDiv, std::vector< std::vector >& intersectEdge) throw(INTERP_KERNEL::Exception); - static void BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const std::vector& b1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, + static void BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, const std::vector< std::vector >& colinear2, const MEDCouplingUMesh *m2, const std::vector& b2, const int *desc2, const int *descIndx2, const std::vector >& intesctEdges2, const std::vector& addCoords, std::vector& cr, std::vector& crI, std::vector& cNb1, std::vector& cNb2); diff --git a/src/MEDLoader/MEDFileMesh.cxx b/src/MEDLoader/MEDFileMesh.cxx index f0b4c0769..cd2284a5a 100644 --- a/src/MEDLoader/MEDFileMesh.cxx +++ b/src/MEDLoader/MEDFileMesh.cxx @@ -1166,6 +1166,14 @@ int MEDFileUMesh::getMeshDimension() const throw(INTERP_KERNEL::Exception) throw INTERP_KERNEL::Exception("MEDFileUMesh::getMeshDimension : impossible to find a mesh dimension !"); } +int MEDFileUMesh::getSpaceDimension() const throw(INTERP_KERNEL::Exception) +{ + const DataArrayDouble *coo=_coords; + if(!coo) + throw INTERP_KERNEL::Exception(" MEDFileUMesh::getSpaceDimension : no coords set !"); + return coo->getNumberOfComponents(); +} + std::string MEDFileUMesh::simpleRepr() const { std::ostringstream oss; diff --git a/src/MEDLoader/MEDFileMesh.hxx b/src/MEDLoader/MEDFileMesh.hxx index 013fcb2b7..8aa496ecf 100644 --- a/src/MEDLoader/MEDFileMesh.hxx +++ b/src/MEDLoader/MEDFileMesh.hxx @@ -142,6 +142,7 @@ namespace ParaMEDMEM // void write(const char *fileName, int mode) const throw(INTERP_KERNEL::Exception); int getMeshDimension() const throw(INTERP_KERNEL::Exception); + int getSpaceDimension() const throw(INTERP_KERNEL::Exception); std::string simpleRepr() const; std::string advancedRepr() const; int getSizeAtLevel(int meshDimRelToMaxExt) const throw(INTERP_KERNEL::Exception); diff --git a/src/MEDLoader/Swig/MEDLoader.i b/src/MEDLoader/Swig/MEDLoader.i index ec3ab77a0..b10981ab5 100644 --- a/src/MEDLoader/Swig/MEDLoader.i +++ b/src/MEDLoader/Swig/MEDLoader.i @@ -458,6 +458,7 @@ namespace ParaMEDMEM static MEDFileUMesh *New(const char *fileName) throw(INTERP_KERNEL::Exception); static MEDFileUMesh *New(); ~MEDFileUMesh(); + int getSpaceDimension() const throw(INTERP_KERNEL::Exception); // std::vector getGrpNonEmptyLevels(const char *grp) const throw(INTERP_KERNEL::Exception); std::vector getGrpNonEmptyLevelsExt(const char *grp) const throw(INTERP_KERNEL::Exception); -- 2.39.2