From 8d332f198fef163966d75b9b0f205680dad2241d Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 14 Oct 2011 10:33:49 +0000 Subject: [PATCH] Partitionnizer of meshes. --- .../Bases/NormalizedUnstructuredMesh.hxx | 3 +- src/INTERP_KERNEL/CellModel.cxx | 107 ++++- src/INTERP_KERNEL/CellModel.hxx | 1 + .../InterpKernelGeo2DComposedEdge.cxx | 34 +- .../InterpKernelGeo2DComposedEdge.hxx | 2 + .../Geometric2D/InterpKernelGeo2DEdge.cxx | 28 +- .../Geometric2D/InterpKernelGeo2DEdge.hxx | 24 +- .../Geometric2D/InterpKernelGeo2DEdge.txx | 2 +- .../InterpKernelGeo2DEdgeArcCircle.cxx | 41 ++ .../InterpKernelGeo2DEdgeArcCircle.hxx | 5 + .../Geometric2D/InterpKernelGeo2DEdgeLin.cxx | 81 ++++ .../Geometric2D/InterpKernelGeo2DEdgeLin.hxx | 5 + .../InterpKernelGeo2DElementaryEdge.cxx | 29 ++ .../InterpKernelGeo2DElementaryEdge.hxx | 6 + .../Geometric2D/InterpKernelGeo2DNode.cxx | 51 +++ .../Geometric2D/InterpKernelGeo2DNode.hxx | 6 + .../InterpKernelGeo2DQuadraticPolygon.cxx | 145 +++++- .../InterpKernelGeo2DQuadraticPolygon.hxx | 6 + .../InterpKernelCellSimplify.cxx | 35 +- .../InterpKernelCellSimplify.hxx | 2 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 427 ++++++++++++++++-- src/MEDCoupling/MEDCouplingUMesh.hxx | 25 +- src/MEDCoupling_Swig/MEDCoupling.i | 34 ++ 23 files changed, 1018 insertions(+), 81 deletions(-) diff --git a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx index 1a8101da6..a5520aeb3 100644 --- a/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx +++ b/src/INTERP_KERNEL/Bases/NormalizedUnstructuredMesh.hxx @@ -39,6 +39,7 @@ namespace INTERP_KERNEL NORM_POLYGON = 5, NORM_TRI6 = 6, NORM_QUAD8 = 8, + NORM_QPOLYG = 32, // NORM_TETRA4 = 14, NORM_PYRA5 = 15, @@ -51,7 +52,7 @@ namespace INTERP_KERNEL NORM_HEXA20 = 30, NORM_POLYHED = 31, NORM_ERROR = 40, - NORM_MAXTYPE = 31 // = NORM_POLYHED + NORM_MAXTYPE = 32 } NormalizedCellType; class GenericMesh diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx index 4def4b4a6..f9199b56f 100644 --- a/src/INTERP_KERNEL/CellModel.cxx +++ b/src/INTERP_KERNEL/CellModel.cxx @@ -23,6 +23,7 @@ #include #include +#include #include namespace INTERP_KERNEL @@ -33,7 +34,7 @@ namespace INTERP_KERNEL "NORM_PYRA5", "NORM_PENTA6", "", "NORM_HEXA8", "",//15->19 "NORM_TETRA10", "", "NORM_HEXGP12", "NORM_PYRA13", "",//20->24 "NORM_PENTA15", "", "", "", "",//25->29 - "NORM_HEXA20", "NORM_POLYHED", "", "", "",//30->34 + "NORM_HEXA20", "NORM_POLYHED", "NORM_QPOLYG", "", "",//30->34 "", "", "", "", "",//35->39 "NORM_ERROR"}; @@ -96,6 +97,7 @@ namespace INTERP_KERNEL _map_of_unique_instance.insert(std::make_pair(NORM_HEXA20,CellModel(NORM_HEXA20))); _map_of_unique_instance.insert(std::make_pair(NORM_POLYGON,CellModel(NORM_POLYGON))); _map_of_unique_instance.insert(std::make_pair(NORM_POLYHED,CellModel(NORM_POLYHED))); + _map_of_unique_instance.insert(std::make_pair(NORM_QPOLYG,CellModel(NORM_QPOLYG))); } CellModel::CellModel(NormalizedCellType type):_type(type) @@ -274,6 +276,11 @@ namespace INTERP_KERNEL _nb_of_pts=0; _nb_of_sons=0; _dim=3; _dyn=true; _is_simplex=false; } break; + case NORM_QPOLYG: + { + _nb_of_pts=0; _nb_of_sons=0; _dim=2; _dyn=true; _is_simplex=false; + } + break; case NORM_ERROR: { _nb_of_pts=std::numeric_limits::max(); _nb_of_sons=std::numeric_limits::max(); _dim=std::numeric_limits::max(); @@ -289,8 +296,13 @@ namespace INTERP_KERNEL { if(!isDynamic()) return getNumberOfSons(); - if(_dim==2)//polygon - return lgth; + if(_dim==2) + { + if(_type==NORM_POLYGON) + return lgth; + else + return lgth/2; + } else return std::count(conn,conn+lgth,-1)+1; } @@ -302,8 +314,13 @@ namespace INTERP_KERNEL { if(!isDynamic()) return getSonType(sonId); - if(_dim==2)//polygon - return NORM_SEG2; + if(_dim==2) + { + if(_type==NORM_POLYGON) + return NORM_SEG2; + else + return NORM_SEG3; + } //polyedron return NORM_POLYGON; } @@ -329,9 +346,19 @@ namespace INTERP_KERNEL { if(_dim==2)//polygon { - sonNodalConn[0]=nodalConn[sonId]; - sonNodalConn[1]=nodalConn[(sonId+1)%lgth]; - return 2; + if(_type==NORM_POLYGON) + { + sonNodalConn[0]=nodalConn[sonId]; + sonNodalConn[1]=nodalConn[(sonId+1)%lgth]; + return 2; + } + else + { + sonNodalConn[0]=nodalConn[sonId]; + sonNodalConn[1]=nodalConn[(sonId+1)%lgth]; + sonNodalConn[2]=nodalConn[sonId+lgth]; + return 3; + } } else {//polyedron @@ -361,7 +388,10 @@ namespace INTERP_KERNEL if(_dim==2)//polygon { - return 2; + if(_type==NORM_POLYGON) + return 2; + else + return 3; } else {//polyedron @@ -376,4 +406,63 @@ namespace INTERP_KERNEL } } + /*! + * This method retrieves if cell1 represented by 'conn1' and cell2 represented by 'conn2' + * are equivalent by a permutation or not. This method expects to work on 1D or 2D (only mesh dimension where it is possible to have a spaceDim) strictly higher than meshDim. + * If not an exception will be thrown. + * @return True if two cells have same orientation, false if not. + */ + bool CellModel::getOrientationStatus(unsigned lgth, const int *conn1, const int *conn2) const + { + if(_dim!=1 && _dim!=2) + throw INTERP_KERNEL::Exception("CellModel::getOrientationStatus : invalid dimension ! Must be 1 or 2 !"); + if(!_quadratic) + { + std::vector tmp(2*lgth); + std::vector::iterator it=std::copy(conn1,conn1+lgth,tmp.begin()); + std::copy(conn1,conn1+lgth,it); + it=std::find_first_of(tmp.begin(),tmp.end(),conn2,conn2+lgth); + return it!=tmp.end(); + } + else + { + if(_dim!=1) + { + std::vector tmp(lgth); + std::vector::iterator it=std::copy(conn1,conn1+lgth/2,tmp.begin()); + std::copy(conn1,conn1+lgth/2,it); + it=std::find_first_of(tmp.begin(),tmp.end(),conn2,conn2+lgth/2); + int d=std::distance(tmp.begin(),it); + if(it==tmp.end()) + return false; + it=std::copy(conn1+lgth/2,conn1+lgth,tmp.begin()); + std::copy(conn1+lgth/2,conn1+lgth,it); + it=std::find_first_of(tmp.begin(),tmp.end(),conn2,conn2+lgth); + if(it==tmp.end()) + return false; + int d2=std::distance(tmp.begin(),it); + return d==d2; + } + else + { + int p=(lgth+1)/2; + std::vector tmp(2*p); + std::vector::iterator it=std::copy(conn1,conn1+p,tmp.begin()); + std::copy(conn1,conn1+p,it); + it=std::find_first_of(tmp.begin(),tmp.end(),conn2,conn2+p); + int d=std::distance(tmp.begin(),it); + if(it==tmp.end()) + return false; + tmp.resize(2*p-2); + it=std::copy(conn1+p,conn1+lgth,tmp.begin()); + std::copy(conn1+p,conn1+lgth,it); + it=std::find_first_of(tmp.begin(),tmp.end(),conn2+p,conn2+lgth); + if(it==tmp.end()) + return false; + int d2=std::distance(tmp.begin(),it); + return d==d2; + } + } + } + } diff --git a/src/INTERP_KERNEL/CellModel.hxx b/src/INTERP_KERNEL/CellModel.hxx index 8c4b7d236..4a0706646 100644 --- a/src/INTERP_KERNEL/CellModel.hxx +++ b/src/INTERP_KERNEL/CellModel.hxx @@ -50,6 +50,7 @@ namespace INTERP_KERNEL INTERPKERNEL_EXPORT bool isSimplex() const { return _is_simplex; } //! sonId is in C format. INTERPKERNEL_EXPORT const unsigned *getNodesConstituentTheSon(unsigned sonId) const { return _sons_con[sonId]; } + INTERPKERNEL_EXPORT bool getOrientationStatus(unsigned lgth, const int *conn1, const int *conn2) const; INTERPKERNEL_EXPORT unsigned getNumberOfNodes() const { return _nb_of_pts; } INTERPKERNEL_EXPORT unsigned getNumberOfSons() const { return _nb_of_sons; } INTERPKERNEL_EXPORT unsigned getNumberOfSons2(const int *conn, int lgth) const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx index 07b6d5ca3..60f1fa528 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.cxx @@ -222,6 +222,18 @@ double ComposedEdge::normalize(ComposedEdge *other, double& xBary, double& yBary return dimChar; } +double ComposedEdge::normalizeExt(ComposedEdge *other, double& xBary, double& yBary) +{ + Bounds b; + b.prepareForAggregation(); + fillBounds(b); + other->fillBounds(b); + double dimChar=b.getCaracteristicDim(); + b.getBarycenter(xBary,yBary); + applyGlobalSimilarity2(other,xBary,yBary,dimChar); + return dimChar; +} + void ComposedEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const { stream.precision(10); @@ -277,6 +289,26 @@ void ComposedEdge::applyGlobalSimilarity(double xBary, double yBary, double dimC (*iter)->applySimilarity(xBary,yBary,dimChar); } +/*! + * Perform Similarity transformation on all elements of this Nodes and Edges on 'this' and 'other'. + * Nodes can be shared between 'this' and 'other'. + */ +void ComposedEdge::applyGlobalSimilarity2(ComposedEdge *other, double xBary, double yBary, double dimChar) +{ + std::set allNodes; + getAllNodes(allNodes); + std::set allNodes2; + other->getAllNodes(allNodes2); + for(std::set::const_iterator it=allNodes2.begin();it!=allNodes2.end();it++) + if(allNodes.find(*it)!=allNodes.end()) + (*it)->declareOn(); + allNodes.insert(allNodes2.begin(),allNodes2.end()); + for(std::set::iterator iter=allNodes.begin();iter!=allNodes.end();iter++) + (*iter)->applySimilarity(xBary,yBary,dimChar); + for(std::list::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++) + (*iter)->applySimilarity(xBary,yBary,dimChar); +} + /*! * This method append to param 'partConsidered' the part of length of subedges IN or ON. * @param partConsidered INOUT param. @@ -360,7 +392,7 @@ bool ComposedEdge::isInOrOut(Node *nodeToTest) const if(val) { Edge *e=val->getPtr(); - std::auto_ptr intersc(Edge::buildIntersectorWith(e1,e)); + std::auto_ptr intersc(Edge::BuildIntersectorWith(e1,e)); bool obviousNoIntersection,areOverlapped; intersc->areOverlappedOrOnlyColinears(0,obviousNoIntersection,areOverlapped); if(obviousNoIntersection) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx index 5e4299c0a..3f0bd9c62 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DComposedEdge.hxx @@ -55,9 +55,11 @@ namespace INTERP_KERNEL void getBarycenter(double *bary) const; void getBarycenterGeneral(double *bary) const; double normalize(ComposedEdge *other, double& xBary, double& yBary); + double normalizeExt(ComposedEdge *other, double& xBary, double& yBary); void fillBounds(Bounds& output) const; void applySimilarity(double xBary, double yBary, double dimChar); void applyGlobalSimilarity(double xBary, double yBary, double dimChar); + void applyGlobalSimilarity2(ComposedEdge *other, double xBary, double yBary, double dimChar); void dispatchPerimeter(double& partConsidered) const; void dispatchPerimeterExcl(double& partConsidered, double& commonPart) const; double dispatchPerimeterAdv(const ComposedEdge& father, std::vector& result) const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx index 509694609..6b06899a4 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.cxx @@ -541,12 +541,12 @@ void Edge::getNormalVector(double *vectOutput) const vectOutput[1]=-tmp; } -Edge *Edge::buildEdgeFrom(Node *start, Node *end) +Edge *Edge::BuildEdgeFrom(Node *start, Node *end) { return new EdgeLin(start,end); } -Edge *Edge::buildFromXfigLine(std::istream& str) +Edge *Edge::BuildFromXfigLine(std::istream& str) { unsigned char type; str >> type; @@ -577,13 +577,13 @@ bool Edge::intersectWith(const Edge *other, MergePoints& commonNode, return false; delete merge; merge=0; - EdgeIntersector *intersector=buildIntersectorWith(this,other); - ret=intersect(this,other,intersector,merge,commonNode,outVal1,outVal2); + EdgeIntersector *intersector=BuildIntersectorWith(this,other); + ret=Intersect(this,other,intersector,merge,commonNode,outVal1,outVal2); delete intersector; return ret; } -bool Edge::intersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode, +bool Edge::IntersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2) { bool rev=intersector->haveTheySameDirection(); @@ -591,14 +591,14 @@ bool Edge::intersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector * Node *f2End=f2->getNode(rev?END:START); TypeOfLocInEdge place1, place2; intersector->getPlacements(f2Start,f2End,place1,place2,commonNode); - int codeForIntersectionCase=combineCodes(place1,place2); - return splitOverlappedEdges(f1,f2,f2Start,f2End,rev,codeForIntersectionCase,outValForF1,outValForF2); + int codeForIntersectionCase=CombineCodes(place1,place2); + return SplitOverlappedEdges(f1,f2,f2Start,f2End,rev,codeForIntersectionCase,outValForF1,outValForF2); } /*! * Perform 1D linear interpolation. Warning distrib1 and distrib2 are expected to be in ascending mode. */ -void Edge::interpolate1DLin(const std::vector& distrib1, const std::vector& distrib2, std::map >& result) +void Edge::Interpolate1DLin(const std::vector& distrib1, const std::vector& distrib2, std::map >& result) { int nbOfV1=distrib1.size()-1; int nbOfV2=distrib2.size()-1; @@ -622,7 +622,7 @@ void Edge::interpolate1DLin(const std::vector& distrib1, const std::vect SegSegIntersector inters(*e1,*e2); bool b1,b2; inters.areOverlappedOrOnlyColinears(0,b1,b2); - if(intersectOverlapped(e1,e2,&inters,commonNode,*f1,*f2)) + if(IntersectOverlapped(e1,e2,&inters,commonNode,*f1,*f2)) { result[i][j]=f1->getCommonLengthWith(*f2)/e1->getCurveLength(); } @@ -635,7 +635,7 @@ void Edge::interpolate1DLin(const std::vector& distrib1, const std::vect n1->decrRef(); n2->decrRef(); n3->decrRef(); n4->decrRef(); } -EdgeIntersector *Edge::buildIntersectorWith(const Edge *e1, const Edge *e2) +EdgeIntersector *Edge::BuildIntersectorWith(const Edge *e1, const Edge *e2) { EdgeIntersector *ret=0; const EdgeLin *tmp1=0; @@ -671,14 +671,14 @@ void Edge::applySimilarity(double xBary, double yBary, double dimChar) _bounds.applySimilarity(xBary,yBary,dimChar); } -bool Edge::intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, +bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2) { bool obviousNoIntersection; bool areOverlapped; intersector->areOverlappedOrOnlyColinears(whereToFind,obviousNoIntersection,areOverlapped); if(areOverlapped) - return intersectOverlapped(f1,f2,intersector,commonNode,outValForF1,outValForF2); + return IntersectOverlapped(f1,f2,intersector,commonNode,outValForF1,outValForF2); if(obviousNoIntersection) return false; std::vector newNodes; @@ -712,7 +712,7 @@ bool Edge::intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersecto return false; } -int Edge::combineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2) +int Edge::CombineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2) { int ret=(int)code1; ret*=OFFSET_FOR_TYPEOFLOCINEDGE; @@ -729,7 +729,7 @@ int Edge::combineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2) * @param direction is param that specifies if e2 and e1 have same directions (true) or opposed (false). * @param code is the code returned by method Edge::combineCodes. */ -bool Edge::splitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *nE, bool direction, int code, +bool Edge::SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *nE, bool direction, int code, ComposedEdge& outVal1, ComposedEdge& outVal2) { Edge *tmp; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx index 28207b13d..0c539f596 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.hxx @@ -219,11 +219,11 @@ namespace INTERP_KERNEL bool changeEndNodeWithAndKeepTrack(Node *otherEndNode, std::vector& track) const; void addSubEdgeInVector(Node *start, Node *end, ComposedEdge& vec) const; void getNormalVector(double *vectOutput) const; - static EdgeIntersector *buildIntersectorWith(const Edge *e1, const Edge *e2); - static Edge *buildFromXfigLine(std::istream& str); - static Edge *buildEdgeFrom(Node *start, Node *end); + static EdgeIntersector *BuildIntersectorWith(const Edge *e1, const Edge *e2); + static Edge *BuildFromXfigLine(std::istream& str); + static Edge *BuildEdgeFrom(Node *start, Node *end); template - static Edge *buildEdgeFrom(Node *start, Node *middle, Node *end); + static Edge *BuildEdgeFrom(Node *start, Node *middle, Node *end); virtual void update(Node *m) = 0; //! returns area between this and axe Ox delimited along Ox by _start and _end. virtual double getAreaOfZone() const = 0; @@ -250,20 +250,26 @@ namespace INTERP_KERNEL const EdgeArcCircle * &arcSeg) const = 0; bool intersectWith(const Edge *other, MergePoints& commonNode, ComposedEdge& outVal1, ComposedEdge& outVal2) const; - static bool intersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode, + static bool IntersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2); - static void interpolate1DLin(const std::vector& distrib1, const std::vector& distrib2, + static void Interpolate1DLin(const std::vector& distrib1, const std::vector& distrib2, std::map >& result); virtual void dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const = 0; bool isEqual(const Edge& other) const; + public: + virtual void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + 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; protected: Edge():_cnt(1),_loc(FULL_UNKNOWN),_start(0),_end(0) { } virtual ~Edge(); - static int combineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2); - static bool intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, + static int CombineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2); + static bool Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2); //! 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, + 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: diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.txx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.txx index c6d3c20b5..bddb26aee 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.txx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.txx @@ -22,7 +22,7 @@ #include "InterpKernelGeo2DEdgeArcCircle.hxx" template -INTERP_KERNEL::Edge *INTERP_KERNEL::Edge::buildEdgeFrom(Node *start, Node *middle, Node *end) +INTERP_KERNEL::Edge *INTERP_KERNEL::Edge::BuildEdgeFrom(Node *start, Node *middle, Node *end) { return new INTERP_KERNEL::EdgeArcCircle(start,middle,end); } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index c33455624..fc978933b 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -21,6 +21,7 @@ #include "InterpKernelGeo2DEdgeLin.hxx" #include "InterpKernelException.hxx" #include "InterpKernelGeo2DNode.hxx" +#include "NormalizedUnstructuredMesh.hxx" #include #include @@ -693,3 +694,43 @@ void EdgeArcCircle::updateBounds() if(isIn2Pi(_angle0,_angle,M_PI)) _bounds[0]=_center[0]-_radius; } + +void EdgeArcCircle::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const +{ + int tmp[2]; + _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); + _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); + if(direction) + { + edgesThis.push_back(tmp[0]); + edgesThis.push_back(tmp[1]); + } + else + { + edgesThis.push_back(tmp[1]); + edgesThis.push_back(tmp[0]); + } +} + +void EdgeArcCircle::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 +{ + _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); + _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); +} + +void EdgeArcCircle::sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis) +{ + Bounds b; + b.prepareForAggregation(); + b.aggregate(getBounds()); + double xBary,yBary; + double dimChar=b.getCaracteristicDim(); + b.getBarycenter(xBary,yBary); + applySimilarity(xBary,yBary,dimChar); + for(std::vector::const_iterator iter=addNodes.begin();iter!=addNodes.end();iter++) + (*iter)->applySimilarity(xBary,yBary,dimChar); + _start->applySimilarity(xBary,yBary,dimChar); + _end->applySimilarity(xBary,yBary,dimChar); +} diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index 8cf60482d..5d2e80f5f 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -111,6 +111,11 @@ namespace INTERP_KERNEL protected: void updateBounds(); Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const; + void fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& 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; + void sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis); protected: //!Value between -2Pi and 2Pi double _angle; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx index 55938f25c..e31aa3cc5 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx @@ -20,6 +20,7 @@ #include "InterpKernelGeo2DEdgeLin.hxx" #include "InterpKernelGeo2DNode.hxx" #include "InterpKernelException.hxx" +#include "NormalizedUnstructuredMesh.hxx" using namespace INTERP_KERNEL; @@ -292,3 +293,83 @@ double EdgeLin::getCharactValueEng(const double *node) const double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1]; return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y); } + +void EdgeLin::fillGlobalInfoAbs(bool direction, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const +{ + int tmp[2]; + _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp); + _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1); + if(direction) + { + edgesThis.push_back(tmp[0]); + edgesThis.push_back(tmp[1]); + } + else + { + edgesThis.push_back(tmp[1]); + edgesThis.push_back(tmp[0]); + } +} + +void EdgeLin::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 +{ + _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); + _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther); +} + +inline bool eqpair(const std::pair& p1, const std::pair& p2) +{ + return fabs(p1.first-p2.first)& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis) +{ + Bounds b; + b.prepareForAggregation(); + b.aggregate(getBounds()); + double xBary,yBary; + double dimChar=b.getCaracteristicDim(); + b.getBarycenter(xBary,yBary); + for(std::vector::const_iterator iter=addNodes.begin();iter!=addNodes.end();iter++) + (*iter)->applySimilarity(xBary,yBary,dimChar); + _start->applySimilarity(xBary,yBary,dimChar); + _end->applySimilarity(xBary,yBary,dimChar); + std::size_t sz=addNodes.size(); + std::vector< std::pair > an2(sz); + for(std::size_t i=0;i(getCharactValue(*addNodes[i]),addNodes[i]); + std::sort(an2.begin(),an2.end()); + int startId=(*mapp1.find(_start)).second; + int endId=(*mapp1.find(_end)).second; + std::vector tmpp; + std::vector< std::pair >::const_iterator itend=std::unique(an2.begin(),an2.end(),eqpair); + for(std::vector< std::pair >::const_iterator it=an2.begin();it!=itend;it++) + { + int idd=(*mapp2.find((*it).second)).second; + if((*it).first1-QUADRATIC_PLANAR::_precision) + { + endId=idd; + continue; + } + tmpp.push_back(idd); + } + std::vector tmpp2(tmpp.size()+2); + tmpp2[0]=startId; + std::copy(tmpp.begin(),tmpp.end(),tmpp2.begin()+1); + tmpp2[tmpp.size()+1]=endId; + std::vector::iterator itt=std::unique(tmpp2.begin(),tmpp2.end()); + tmpp2.resize(std::distance(tmpp2.begin(),itt)); + int nbOfEdges=tmpp2.size()-1; + for(int i=0;i& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& 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; + void sortIdsAbs(const std::vector& addNodes, const std::map& mapp1, const std::map& mapp2, std::vector& edgesThis); }; } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx index 4369f2486..ff5a98a4f 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.cxx @@ -200,3 +200,32 @@ bool ElementaryEdge::isEqual(const ElementaryEdge& other) const { return _ptr->isEqual(*other._ptr); } + +/*! + * Called by QuadraticPolygon::splitAbs method. + */ +void ElementaryEdge::fillGlobalInfoAbs(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& edgesThis, std::vector& addCoo, std::map mapAddCoo) const +{ + _ptr->fillGlobalInfoAbs(_direction,mapThis,mapOther,offset1,offset2,fact,baryX,baryY,edgesThis,addCoo,mapAddCoo); +} + +/*! + * Called by QuadraticPolygon::splitAbs method. Close to ElementaryEdge::fillGlobalInfoAbs method expect that here edgesOther (that replace edgesThis) is here an in/out parameter that only contains nodes + * unsorted because the "other" mesh is not subdivided yet. + */ +void ElementaryEdge::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 +{ + _ptr->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,edgesOther,addCoo,mapAddCoo); +} + +/*! + * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the + * orientation of edge. Called by QuadraticPolygon::buildFromCrudeDataArray. + */ +ElementaryEdge *ElementaryEdge::BuildEdgeFromCrudeDataArray(bool isQuad, bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end) +{ + Edge *ptr=Edge::BuildEdgeFrom(start,end); + return new ElementaryEdge(ptr,direction); +} diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx index e70dfd0ba..5a23a6d96 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DElementaryEdge.hxx @@ -65,6 +65,12 @@ namespace INTERP_KERNEL bool getDirection() const { return _direction; } bool intresincEqCoarse(const Edge *other) const; bool isEqual(const ElementaryEdge& other) const; + public: + void fillGlobalInfoAbs(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& 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); private: bool _direction; Edge *_ptr; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx index 9ac557b0e..7e336ade2 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx @@ -130,3 +130,54 @@ void Node::applySimilarity(double xBary, double yBary, double dimChar) _coords[0]=(_coords[0]-xBary)/dimChar; _coords[1]=(_coords[1]-yBary)/dimChar; } + +/*! + * Called by QuadraticPolygon::splitAbs method. + */ +void Node::fillGlobalInfoAbs(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& addCoo, std::map mapAddCoo, int *nodeId) const +{ + std::map::const_iterator it=mapThis.find(const_cast(this)); + if(it!=mapThis.end()) + { + *nodeId=(*it).second; + return; + } + it=mapOther.find(const_cast(this)); + if(it!=mapOther.end()) + { + *nodeId=(*it).second+offset1; + return; + } + it=mapAddCoo.find(const_cast(this)); + if(it!=mapAddCoo.end()) + { + *nodeId=(*it).second; + return; + } + int id=addCoo.size()/2; + addCoo.push_back(fact*_coords[0]+baryX); + addCoo.push_back(fact*_coords[0]+baryX); + mapAddCoo[const_cast(this)]=offset2+id; +} + +/*! + * Called by QuadraticPolygon::splitAbs method. + */ +void Node::fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& addCoo, std::map mapAddCoo, std::vector& pointsOther) const +{ + int tmp; + std::size_t sz1=addCoo.size(); + fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,&tmp); + if(sz1!=addCoo.size()) + { + pointsOther.push_back(tmp); + return ; + } + std::vector::const_iterator it=std::find(pointsOther.begin(),pointsOther.end(),tmp); + if(it!=pointsOther.end()) + return ; + if(tmp #include #include #include @@ -83,6 +84,11 @@ namespace INTERP_KERNEL static bool areDoubleEqualsWP(double a, double b, double k) { return fabs(a-b) < k*QUADRATIC_PLANAR::_precision; } static double distanceBtw2Pt(const double *a, const double *b) { return sqrt((a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1])); } static double distanceBtw2PtSq(const double *a, const double *b) { return (a[0]-b[0])*(a[0]-b[0])+(a[1]-b[1])*(a[1]-b[1]); } + // + void fillGlobalInfoAbs(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& addCoo, std::map mapAddCoo, int *nodeId) const; + void fillGlobalInfoAbs2(const std::map& mapThis, const std::map& mapOther, int offset1, int offset2, double fact, double baryX, double baryY, + std::vector& addCoo, std::map mapAddCoo, std::vector& pointsOther) const; protected: ~Node(); protected: diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index d66bdae4e..845e9fc5d 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -49,7 +49,7 @@ QuadraticPolygon::QuadraticPolygon(const char *file) while(strcmp(currentLine,"1200 2")!=0); do { - Edge *newEdge=Edge::buildFromXfigLine(stream); + Edge *newEdge=Edge::BuildFromXfigLine(stream); if(!empty()) newEdge->changeStartNodeWith(back()->getEndNode()); pushBack(newEdge); @@ -207,6 +207,149 @@ double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other) return ret*fact*fact; } +/*! + * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondance between nodes contained in 'this' and node ids in a global mesh. + * In the same way, 'mapOther' gives the correspondance between nodes contained in 'other' and node ids in a global mesh from wich 'other' is extracted. + * This method has 1 out paramater : 'edgesThis', After the call of this method contains nodal connectivity (including type) of 'this' into globlal "this mesh". + * This method has 2 in/out parameters : 'subDivOther' and 'addCoo'.'otherEdgeIds' is useful to put values in 'edgesThis', 'subDivOther' and 'addCoo'. + * Size of 'otherEdgeIds' has to be equal to number of ElementaryEdges in 'other'. No check of that will be done. + * @param offset1 is the number of nodes contained in global mesh from which 'this' is extracted. + * @param offset2 is the sum of nodes contained in global mesh from which 'this' is extracted and 'other' is extracted. + * @otherEdgeIds is a vector with the same size than other before calling this method. It gives in the same order the cell id in global other mesh. + */ +void QuadraticPolygon::splitAbs(QuadraticPolygon& other, const std::map& mapThis, const std::map& mapOther, int offset1, int offset2 , const std::vector& otherEdgeIds, + std::vector& edgesThis, std::vector< std::vector >& subDivOther, std::vector& addCoo) +{ + double xBaryBB, yBaryBB; + double fact=normalize(&other, xBaryBB, yBaryBB); + // + IteratorOnComposedEdge it1(this),it3(&other); + MergePoints merge; + ComposedEdge *c1=new ComposedEdge; + ComposedEdge *c2=new ComposedEdge; + int i=0; + std::map mapAddCoo; + for(it3.first();!it3.finished();it3.next(),i++) + { + QuadraticPolygon otherTmp; + ElementaryEdge* curE3=it3.current(); + otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef(); + IteratorOnComposedEdge it2(&otherTmp); + for(it2.first();it2.finished();it2.next()) + { + ElementaryEdge* curE2=it2.current(); + if(!curE2->isThereStartPoint()) + it1.first(); + else + it1=curE2->getIterator(); + for(;!it1.finished();) + { + ElementaryEdge* curE1=it1.current(); + merge.clear(); + if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2)) + { + if(!curE1->getDirection()) c1->reverse(); + if(!curE2->getDirection()) c2->reverse(); + updateNeighbours(merge,it1,it2,c1,c2); + //Substitution of simple edge by sub-edges. + delete curE1; // <-- destroying simple edge coming from pol1 + delete curE2; // <-- destroying simple edge coming from pol2 + it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next. + it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next. + curE2=it2.current(); + // + it1.assignMySelfToAllElems(c2);//To avoid that others + SoftDelete(c1); + SoftDelete(c2); + c1=new ComposedEdge; + c2=new ComposedEdge; + } + else + { + updateNeighbours(merge,it1,it2,curE1,curE2); + it1.next(); + } + } + } + if(otherTmp._sub_edges.size()>1) + { + for(std::list::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++) + (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo); + } + } + Delete(c1); + Delete(c2); + // + for(std::list::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++) + (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo); + // +} + +/*! + * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the + * orientation of edge. + */ +void QuadraticPolygon::buildFromCrudeDataArray(const std::map& mapp, bool isQuad, const int *descBg, const int *descEnd, const std::vector >& intersectEdges) +{ + int nbOfSeg=std::distance(descBg,descEnd); + for(int i=0;i0; + int edgeId=abs(descBg[i])-1; + const std::vector& subEdge=intersectEdges[edgeId]; + int nbOfSubEdges=subEdge.size()-1; + for(int j=0;j& mapp, std::vector& conn, std::vector& connI) +{ + int nbOfNodesInPg=0,i=0; + for(std::list::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,i++) + { + Node *tmp=0; + if(i==0) + { + tmp=(*it)->getStartNode(); + std::map::const_iterator it=mapp.find(tmp); + conn.push_back((*it).second); + nbOfNodesInPg++; + } + tmp=(*it)->getEndNode(); + std::map::const_iterator it=mapp.find(tmp); + conn.push_back((*it).second); + nbOfNodesInPg++; + } + connI.push_back(connI.back()+nbOfNodesInPg); +} + +/*! + * This method make the hypothesis that 'this' and 'other' are splited at the minimum into edges that are fully IN, OUT or ON. + * This method returns newly created polygons in 'conn' and 'connI' and the corresponding ids ('idThis','idOther') are stored respectively into 'nbThis' and 'nbOther'. + */ +void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, const std::map& mapp, int idThis, int idOther, std::vector& conn, std::vector& connI, std::vector& nbThis, std::vector& nbOther) +{ + double xBaryBB, yBaryBB; + normalizeExt(&other, xBaryBB, yBaryBB); + //Locate 'this' relative to 'other' + other.performLocatingOperation(*this); + dumpInXfigFileWithOther(other,"tony.fig"); + std::vector res=other.buildIntersectionPolygons(*this,other); + for(std::vector::iterator it=res.begin();it!=res.end();it++) + { + (*it)->appendCrudeData(mapp,conn,connI); + nbThis.push_back(idThis); + nbOther.push_back(idOther); + delete *it; + } +} + /*! * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method. * 'other' is a QuadraticPolygon of \b non closed edges. diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index 5b9e2bc6e..dcfe0ecb2 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -54,6 +54,12 @@ namespace INTERP_KERNEL double intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear); //! Before intersecting as intersectWith a normalization is done. 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, 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 appendCrudeData(const std::map& mapp, std::vector& conn, std::vector& connI); + void buildPartitionsAbs(QuadraticPolygon& other, const std::map& mapp, int idThis, int idOther, std::vector& conn, std::vector& connI, std::vector& nb1, std::vector& nb2); + double intersectWith(const QuadraticPolygon& other) const; double intersectWith(const QuadraticPolygon& other, double* barycenter) const; std::vector intersectMySelfWith(const QuadraticPolygon& other) const; diff --git a/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx b/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx index 5ebce3c4c..d8d5e6ae9 100644 --- a/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx +++ b/src/INTERP_KERNEL/InterpKernelCellSimplify.cxx @@ -60,7 +60,7 @@ INTERP_KERNEL::NormalizedCellType CellSimplify::simplifyDegeneratedCell(INTERP_K for(int i=1;i #include @@ -506,11 +510,57 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA } } +/// @cond INTERNAL + +int MEDCouplingFastNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2) +{ + return id; +} + +int MEDCouplingOrientationSensitiveNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2) +{ + if(!compute) + return id+1; + else + { + if(cm.getOrientationStatus(nb,conn1,conn2)) + return id+1; + else + return -(id+1); + } +} + +/// @endcond + /*! * \b WARNING this method do the assumption that connectivity lies on the coordinates set. * For speed reasons no check of this will be done. */ MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception) +{ + return buildDescendingConnectivityGen(desc,descIndx,revDesc,revDescIndx,MEDCouplingFastNbrer); +} + +/*! + * WARNING this method do the assumption that connectivity lies on the coordinates set. + * For speed reasons no check of this will be done. + * This method differs from MEDCouplingUMesh::buildDescendingConnectivity method in that 'desc' is in different format. + * This method is more precise because it returns in descending connectivity giving the direction. If value is positive the n-1 dim element is taken in the same direction, + * if it is in the opposite direction it is retrieved negative. So the problem is for elemt #0 in C convention. That's why this method is the only one that retrieves + * an array in relative "FORTRAN" mode. + */ +MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception) +{ + return buildDescendingConnectivityGen(desc,descIndx,revDesc,revDescIndx,MEDCouplingOrientationSensitiveNbrer); +} + +/// @cond INTERNAL + +/*! + * \b WARNING this method do the assumption that connectivity lies on the coordinates set. + * For speed reasons no check of this will be done. + */ +MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const throw(INTERP_KERNEL::Exception) { checkFullyDefined(); int nbOfCells=getNumberOfCells(); @@ -568,13 +618,13 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivity(DataArrayInt *de revNodalB[tmp[k]].push_back(cellDM1Id); revDescMeshConnB.resize(cellDM1Id+1); revDescMeshConnB.back().push_back(eltId); - descMeshConnB[eltId].push_back(cellDM1Id); + descMeshConnB[eltId].push_back(nbrer(cellDM1Id,0,cms,false,0,0)); } else { int DM1cellId=shareableCellsL.front(); revDescMeshConnB[DM1cellId].push_back(eltId); - descMeshConnB[eltId].push_back(DM1cellId); + descMeshConnB[eltId].push_back(nbrer(DM1cellId,nbOfNodesSon,cms,true,tmp,&meshDM1Conn[meshDM1ConnIndex[DM1cellId]])); } } delete [] tmp; @@ -618,6 +668,9 @@ struct MEDCouplingAccVisit int _new_nb_of_nodes; }; +/// @endcond + + /*! * This method convert this into dynamic types without changing geometry. * That is to say if 'this' is a 2D, mesh after the invocation of this method it will contain only polygons. @@ -637,7 +690,13 @@ void MEDCouplingUMesh::convertToPolyTypes(const std::vector& cellIdsToConve const int *connIndex=_nodal_connec_index->getConstPointer(); int *conn=_nodal_connec->getPointer(); for(std::vector::const_iterator iter=cellIdsToConvert.begin();iter!=cellIdsToConvert.end();iter++) - conn[connIndex[*iter]]=INTERP_KERNEL::NORM_POLYGON; + { + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connIndex[*iter]]); + if(!cm.isDynamic()) + conn[connIndex[*iter]]=INTERP_KERNEL::NORM_POLYGON; + else + conn[connIndex[*iter]]=INTERP_KERNEL::NORM_QPOLYG; + } } else { @@ -800,7 +859,7 @@ void MEDCouplingUMesh::unPolyze() { INTERP_KERNEL::AutoPtr tmp=new int[lgthOfCurCell-1]; std::copy(conn+posOfCurCell+1,conn+posOfCurCell+lgthOfCurCell,(int *)tmp); - newType=INTERP_KERNEL::CellSimplify::tryToUnPoly2D(tmp,lgthOfCurCell-1,conn+newPos+1,newLgth); + newType=INTERP_KERNEL::CellSimplify::tryToUnPoly2D(cm.isQuadratic(),tmp,lgthOfCurCell-1,conn+newPos+1,newLgth); } if(cm.getDimension()==3) { @@ -2370,6 +2429,112 @@ namespace ParaMEDMEM INTERP_KERNEL::NormalizedCellType getTypeOfElement(int) const { return (INTERP_KERNEL::NormalizedCellType)0; } // end }; + + INTERP_KERNEL::Edge *MEDCouplingUMeshBuildQPFromEdge(INTERP_KERNEL::NormalizedCellType typ, std::map& mapp2, const int *bg) + { + INTERP_KERNEL::Edge *ret=0; + switch(typ) + { + case INTERP_KERNEL::NORM_SEG2: + { + ret=new INTERP_KERNEL::EdgeLin(mapp2[bg[0]],mapp2[bg[1]]); + break; + } + case INTERP_KERNEL::NORM_SEG3: + { + ret=new INTERP_KERNEL::EdgeArcCircle(mapp2[bg[0]],mapp2[bg[2]],mapp2[bg[1]]); + break; + } + default: + throw INTERP_KERNEL::Exception("MEDCouplingUMeshBuildQPFromEdge : Expecting a mesh with spaceDim==2 and meshDim==1 !"); + } + return ret; + } + + /*! + * This method creates a sub mesh in Geometric2D DS. The sub mesh is composed be the sub set of cells in 'candidates' and the global mesh 'mDesc'. + * The input meth 'mDesc' must be so that mDim==1 et spaceDim==3. + * 'mapp' contains a mapping between local numbering in submesh and the global node numbering in 'mDesc'. + */ + INTERP_KERNEL::QuadraticPolygon *MEDCouplingUMeshBuildQPFromMesh(const MEDCouplingUMesh *mDesc, const std::vector& candidates, std::map& mapp) throw(INTERP_KERNEL::Exception) + { + mapp.clear(); + std::map mapp2; + const double *coo=mDesc->getCoords()->getConstPointer(); + const int *c=mDesc->getNodalConnectivity()->getConstPointer(); + const int *cI=mDesc->getNodalConnectivityIndex()->getConstPointer(); + std::set s; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + s.insert(c+cI[*it]+1,c+cI[(*it)+1]); + for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) + { + INTERP_KERNEL::Node *n=new INTERP_KERNEL::Node(coo[2*(*it2)],coo[2*(*it2)+1]); + mapp[n]=(*it2); mapp2[*it2]=n; + } + INTERP_KERNEL::QuadraticPolygon *ret=new INTERP_KERNEL::QuadraticPolygon; + for(std::vector::const_iterator it=candidates.begin();it!=candidates.end();it++) + { + INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)c[cI[*it]]; + ret->pushBack(MEDCouplingUMeshBuildQPFromEdge(typ,mapp2,c+cI[*it]+1)); + } + for(std::set::const_iterator it2=s.begin();it2!=s.end();it2++) + mapp2[*it2]->decrRef(); + return ret; + } + + INTERP_KERNEL::Node *MEDCouplingUMeshBuildQPNode(int nodeId, const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo) + { + if(nodeId>=offset2) + { + int locId=nodeId-offset2; + return new INTERP_KERNEL::Node(addCoo[2*locId],addCoo[2*locId+1]); + } + if(nodeId>=offset1) + { + int locId=nodeId-offset1; + return new INTERP_KERNEL::Node(coo2[2*locId],coo2[2*locId+1]); + } + return new INTERP_KERNEL::Node(coo1[2*nodeId],coo1[2*nodeId+1]); + } + + void MEDCouplingUMeshBuildQPFromMesh2(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector& addCoo, + bool isQuad1, const int *desc1Bg, const int *desc1End, const std::vector >& intesctEdges1, + bool isQuad2, const int *desc2Bg, const int *desc2End, const std::vector >& intesctEdges2, + /*output*/INTERP_KERNEL::QuadraticPolygon& pol1, INTERP_KERNEL::QuadraticPolygon& pol2, std::map& mapp) + { + std::map mappRev; + for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++) + { + int eltId1=abs(*desc1)-1; + for(std::vector::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++) + { + std::map::const_iterator it=mappRev.find(*it1); + if(it==mappRev.end()) + { + INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo); + mapp[node]=*it1; + mappRev[*it1]=node; + } + } + } + for(const int *desc2=desc2Bg;desc2!=desc2End;desc2++) + { + int eltId2=abs(*desc2)-1; + for(std::vector::const_iterator it2=intesctEdges2[eltId2].begin();it2!=intesctEdges2[eltId2].end();it2++) + { + std::map::const_iterator it=mappRev.find(*it2); + if(it==mappRev.end()) + { + INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it2,coo1,offset1,coo2,offset2,addCoo); + mapp[node]=*it2; + mappRev[*it2]=node; + } + } + } + // + pol1.buildFromCrudeDataArray(mappRev,isQuad1,desc1Bg,desc1End,intesctEdges1); + pol2.buildFromCrudeDataArray(mappRev,isQuad2,desc2Bg,desc2End,intesctEdges2); + } } /// @endcond @@ -2780,7 +2945,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMeshFromThisLowLev(int nbOfNode std::vector newc; for(int j=0;jalloc(newc.size()*nbOf1DCells,1); @@ -3111,9 +3276,10 @@ void MEDCouplingUMesh::are2DCellsNotCorrectlyOriented(const double *vec, bool po for(int i=0;i tmp(connI[i+1]-connI[i]-2); - std::copy(conn+connI[i]+2,conn+connI[i+1],tmp.rbegin()); - std::copy(tmp.begin(),tmp.end(),conn+connI[i]+2); - } + if(!polyOnly || (type==INTERP_KERNEL::NORM_POLYGON || type==INTERP_KERNEL::NORM_QPOLYG)) + { + bool isQuadratic=INTERP_KERNEL::CellModel::GetCellModel(type).isQuadratic(); + if(!IsPolygonWellOriented(isQuadratic,vec,conn+connI[i]+1,conn+connI[i+1],coordsPtr)) + { + isModified=true; + std::vector tmp(connI[i+1]-connI[i]-2); + std::copy(conn+connI[i]+2,conn+connI[i+1],tmp.rbegin()); + std::copy(tmp.begin(),tmp.end(),conn+connI[i]+2); + } + } } if(isModified) _nodal_connec->declareAsNew(); @@ -3254,19 +3423,19 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getEdgeRatioField() const throw(INTERP { case INTERP_KERNEL::NORM_TRI3: { - fillInCompact3DMode(spaceDim,3,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp); *pt=INTERP_KERNEL::triEdgeRatio(tmp); break; } case INTERP_KERNEL::NORM_QUAD4: { - fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::quadEdgeRatio(tmp); break; } case INTERP_KERNEL::NORM_TETRA4: { - fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::tetraEdgeRatio(tmp); break; } @@ -3314,19 +3483,19 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getAspectRatioField() const throw(INTE { case INTERP_KERNEL::NORM_TRI3: { - fillInCompact3DMode(spaceDim,3,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,3,conn+1,coo,tmp); *pt=INTERP_KERNEL::triAspectRatio(tmp); break; } case INTERP_KERNEL::NORM_QUAD4: { - fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::quadAspectRatio(tmp); break; } case INTERP_KERNEL::NORM_TETRA4: { - fillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); + FillInCompact3DMode(spaceDim,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::tetraAspectRatio(tmp); break; } @@ -3374,7 +3543,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getWarpField() const throw(INTERP_KERN { case INTERP_KERNEL::NORM_QUAD4: { - fillInCompact3DMode(3,4,conn+1,coo,tmp); + FillInCompact3DMode(3,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::quadWarp(tmp); break; } @@ -3422,7 +3591,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getSkewField() const throw(INTERP_KERN { case INTERP_KERNEL::NORM_QUAD4: { - fillInCompact3DMode(3,4,conn+1,coo,tmp); + FillInCompact3DMode(3,4,conn+1,coo,tmp); *pt=INTERP_KERNEL::quadSkew(tmp); break; } @@ -3730,8 +3899,8 @@ MEDCouplingUMesh *MEDCouplingUMesh::emulateMEDMEMBDC(const MEDCouplingUMesh *nM1 */ DataArrayInt *MEDCouplingUMesh::sortCellsInMEDFileFrmt() throw(INTERP_KERNEL::Exception) { - static const int N=18; - static const INTERP_KERNEL::NormalizedCellType MEDMEM_ORDER[N] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_POLYHED }; + static const int N=19; + static const INTERP_KERNEL::NormalizedCellType MEDMEM_ORDER[N] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_POLYHED }; checkConnectivityFullyDefined(); MEDCouplingAutoRefCountObjectPtr ret=getRenumArrForConsecutiveCellTypesSpec(MEDMEM_ORDER,MEDMEM_ORDER+N); renumberCells(ret->getConstPointer(),false); @@ -3986,6 +4155,25 @@ MEDCouplingUMesh *MEDCouplingUMesh::keepSpecifiedCells(INTERP_KERNEL::Normalized return ret2; } +/*! + * This method returns a vector of size 'this->getNumberOfCells()'. + * This method retrieves for each cell in 'this' if it is linear (false) or quadratic(true). + */ +std::vector MEDCouplingUMesh::getQuadraticStatus() const throw(INTERP_KERNEL::Exception) +{ + int ncell=getNumberOfCells(); + std::vector ret(ncell); + const int *cI=getNodalConnectivityIndex()->getConstPointer(); + const int *c=getNodalConnectivity()->getConstPointer(); + for(int i=0;i& ret) +void MEDCouplingUMesh::AppendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector& ret) { INTERP_KERNEL::NormalizedCellType flatType=(INTERP_KERNEL::NormalizedCellType)connBg[0]; const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(flatType); @@ -4334,10 +4522,12 @@ void MEDCouplingUMesh::appendExtrudedCell(const int *connBg, const int *connEnd, /*! * This static operates only for coords in 3D. The polygon is specfied by its connectivity nodes in [begin,end). */ -bool MEDCouplingUMesh::IsPolygonWellOriented(const double *vec, const int *begin, const int *end, const double *coords) +bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords) { double v[3]={0.,0.,0.}; int sz=std::distance(begin,end); + if(isQuadratic) + sz/=2; for(int i=0;i\n"; ofs << " \n"; ofs << " \n" << pointData << std::endl; @@ -4482,7 +4672,7 @@ void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData c2->writeVTK(ofs,8,"UInt8","types"); ofs << " \n"; ofs << " \n"; - ofs << " \n"; + ofs << " \n"; } std::string MEDCouplingUMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception) @@ -4490,6 +4680,181 @@ std::string MEDCouplingUMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exc return std::string("UnstructuredGrid"); } +MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) throw(INTERP_KERNEL::Exception) +{ + std::vector< std::vector > intersectEdge1, subDiv2; + MEDCouplingUMesh *m1Desc=0,*m2Desc=0; + DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0; + std::vector addCoo; + INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps; + IntersectDescending2DMeshes(m1,m2,eps,intersectEdge1,subDiv2,m1Desc,desc1,descIndx1,revDesc1,revDescIndx1, + m2Desc,desc2,descIndx2,revDesc2,revDescIndx2,addCoo); + revDesc1->decrRef(); revDescIndx1->decrRef(); revDesc2->decrRef(); revDescIndx2->decrRef(); + std::vector< std::vector > intersectEdge2; + BuildIntersectEdges(m1Desc,m2Desc,addCoo,subDiv2,intersectEdge2); + std::vector b1=m1Desc->getQuadraticStatus(); + std::vector b2=m1Desc->getQuadraticStatus(); + subDiv2.clear(); m1Desc->decrRef(); m2Desc->decrRef(); + std::vector cr,crI; + std::vector cNb1,cNb2; + BuildIntersecting2DCellsFromEdges(eps,m1,b1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,m2,b2,desc2->getConstPointer(),descIndx2->getConstPointer(),intersectEdge2,addCoo, + /* outputs -> */cr,crI,cNb1,cNb2); +#if 0 + //tony + MEDCouplingAutoRefCountObjectPtr addCooDa=DataArrayDouble::New(); + addCooDa->alloc(addCoo.size()/2,2); + std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer()); + std::vector coordss(3); + coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; + MEDCouplingAutoRefCountObjectPtr coo=DataArrayDouble::Aggregate(coordss); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingUMesh::New("Intersect2D",2); + MEDCouplingAutoRefCountObjectPtr conn=DataArrayInt::New(); conn->alloc(cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer()); + MEDCouplingAutoRefCountObjectPtr connI=DataArrayInt::New(); connI->alloc(crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer()); + MEDCouplingAutoRefCountObjectPtr c1=DataArrayInt::New(); c1->alloc(cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer()); cellNb1=c1; + MEDCouplingAutoRefCountObjectPtr c2=DataArrayInt::New(); c2->alloc(cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer()); cellNb2=c2; + ret->setConnectivity(conn,connI,true); + ret->setCoords(coo); + ret->incrRef(); c1->incrRef(); c2->incrRef(); + return ret; +#endif + desc1->decrRef(); descIndx1->decrRef(); desc2->decrRef(); descIndx2->decrRef(); + return 0; +} + +void MEDCouplingUMesh::BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const std::vector& b1, const int *desc1, const int *descIndx1, const std::vector >& intesctEdges1, + 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) +{ + static const int SPACEDIM=2; + std::vector bbox1,bbox2; + const double *coo1=m1->getCoords()->getConstPointer(); + int offset1=m1->getNumberOfNodes(); + const double *coo2=m2->getCoords()->getConstPointer(); + int offset2=offset1+m2->getNumberOfNodes(); + m1->getBoundingBoxForBBTree(bbox1); + m2->getBoundingBoxForBBTree(bbox2); + BBTree myTree(&bbox2[0],0,0,m2->getNumberOfCells(),-eps); + int ncell1=m1->getNumberOfCells(); + crI.push_back(0); + for(int i=0;i candidates2; + myTree.getIntersectingElems(&bbox1[i*2*SPACEDIM],candidates2); + for(std::vector::const_iterator it2=candidates2.begin();it2!=candidates2.end();it2++) + { + INTERP_KERNEL::QuadraticPolygon pol1,pol2; + std::map mapp; + MEDCouplingUMeshBuildQPFromMesh2(coo1,offset1,coo2,offset2,addCoords, + b1[i],desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1, + b2[*it2],desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2, + /* output */pol1,pol2,mapp); + pol1.buildPartitionsAbs(pol2,mapp,i,*it2,cr,crI,cNb1,cNb2); + } + } +} + +/*! + * This method is private and is the first step of Partition of 2D mesh (spaceDim==2 and meshDim==2). + * + */ +void MEDCouplingUMesh::IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& subDiv2, + MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, + MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2, + std::vector& addCoo) throw(INTERP_KERNEL::Exception) +{ + static const int SPACEDIM=2; + desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New(); revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New(); + desc2=DataArrayInt::New(); descIndx2=DataArrayInt::New(); revDesc2=DataArrayInt::New(); revDescIndx2=DataArrayInt::New(); + m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1); + m2Desc=m2->buildDescendingConnectivity2(desc2,descIndx2,revDesc2,revDescIndx2); + std::vector bbox1,bbox2; + m1Desc->getBoundingBoxForBBTree(bbox1); + m2Desc->getBoundingBoxForBBTree(bbox2); + int ncell1=m1Desc->getNumberOfCells(); + int ncell2=m2Desc->getNumberOfCells(); + intersectEdge1.resize(ncell1); + subDiv2.resize(ncell2); + BBTree myTree(&bbox2[0],0,0,m2Desc->getNumberOfCells(),-eps); + std::vector candidates1(1); + int offset1=m1->getNumberOfNodes(); + int offset2=offset1+m2->getNumberOfNodes(); + for(int i=0;i candidates2; + myTree.getIntersectingElems(&bbox1[i*2*SPACEDIM],candidates2); + if(!candidates2.empty()) + { + std::map map1,map2; + INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2Desc,candidates2,map2); + INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1); + pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],subDiv2,addCoo); + delete pol2; + delete pol1; + } + } +} + +/*! + * This method performs the 2nd step of Partition of 2D mesh. + * This method has 4 inputs : + * - a mesh 'm1' with meshDim==1 and a SpaceDim==2 + * - a mesh 'm2' with meshDim==1 and a SpaceDim==2 + * - subDiv of size 'm->getNumberOfCells()' that lists for each seg cell in 'm' the splitting node ids in randomly sorted. + * The aim of this method is to sort the splitting nodes, if any, and to put in 'intersectEdge' output paramter based on edges of mesh 'm2' + * @param m1 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. Only present for its coords in case of 'subDiv' shares some nodes of 'm1' + * @param m2 is expected to be a mesh of meshDimension equal to 1 and spaceDim equal to 2. No check of that is performed by this method. + * @param addCoo input parameter with additionnal nodes linked to intersection of the 2 meshes. + */ +void MEDCouplingUMesh::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) +{ + int offset1=m1->getNumberOfCells(); + int ncell=m2->getNumberOfCells(); + const int *c=m2->getNodalConnectivity()->getConstPointer(); + const int *cI=m2->getNodalConnectivityIndex()->getConstPointer(); + const double *coo=m2->getCoords()->getConstPointer(); + const double *cooBis=m1->getCoords()->getConstPointer(); + int offset2=offset1+ncell; + intersectEdge.resize(ncell); + for(int i=0;i& divs=subDiv[i]; + int nnode=cI[1]-cI[0]-1; + std::map mapp2; + std::map mapp22; + for(int j=0;j::const_iterator it=mapp2.begin();it!=mapp2.end();it++) + (*it).second->decrRef(); + std::vector addNodes(divs.size()); + std::map mapp3; + for(std::size_t j=0;jsortIdsAbs(addNodes,mapp22,mapp3,intersectEdge[i]); + for(std::vector::const_iterator it=addNodes.begin();it!=addNodes.end();it++) + (*it)->decrRef(); + e->decrRef(); + } +} + MEDCouplingUMeshCellIterator::MEDCouplingUMeshCellIterator(MEDCouplingUMesh *mesh):_mesh(mesh),_cell(new MEDCouplingUMeshCell(mesh)), _own_cell(true),_cell_id(-1),_nb_cell(0) { diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index e6e680c2b..64a22bedf 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -24,6 +24,8 @@ #include "MEDCouplingPointSet.hxx" #include "MEDCouplingMemArray.hxx" +#include "CellModel.hxx" + #include namespace ParaMEDMEM @@ -101,6 +103,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool areCellsIncludedIn(const MEDCouplingUMesh *other, int compType, DataArrayInt *& arr) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT MEDCouplingUMesh *buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT MEDCouplingUMesh *buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *mergeNodes(double precision, bool& areNodesMerged, int& newNbOfNodes); MEDCOUPLING_EXPORT DataArrayInt *mergeNodes2(double precision, bool& areNodesMerged, int& newNbOfNodes); MEDCOUPLING_EXPORT void tryToShareSameCoordsPermute(const MEDCouplingPointSet& other, double epsilon) throw(INTERP_KERNEL::Exception); @@ -163,6 +166,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayInt *keepCellIdsByType(INTERP_KERNEL::NormalizedCellType type, const int *begin, const int *end) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT MEDCouplingUMesh *keepSpecifiedCells(INTERP_KERNEL::NormalizedCellType type, const int *idsPerGeoTypeBg, const int *idsPerGeoTypeEnd) const; + MEDCOUPLING_EXPORT std::vector getQuadraticStatus() const throw(INTERP_KERNEL::Exception); // MEDCOUPLING_EXPORT MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const; MEDCOUPLING_EXPORT DataArrayDouble *getBarycenterAndOwner() const; @@ -173,9 +177,10 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const MEDCouplingUMesh *mesh1, const MEDCouplingUMesh *mesh2) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT static MEDCouplingUMesh *MergeUMeshesOnSameCoords(const std::vector& meshes); MEDCOUPLING_EXPORT static MEDCouplingUMesh *FuseUMeshesOnSameCoords(const std::vector& meshes, int compType, std::vector& corr); - MEDCOUPLING_EXPORT static bool IsPolygonWellOriented(const double *vec, const int *begin, const int *end, const double *coords); + MEDCOUPLING_EXPORT static bool IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords); MEDCOUPLING_EXPORT static bool IsPolyhedronWellOriented(const int *begin, const int *end, const double *coords); MEDCOUPLING_EXPORT static void TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) throw(INTERP_KERNEL::Exception); private: MEDCouplingUMesh(); MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy); @@ -200,8 +205,22 @@ namespace ParaMEDMEM template void getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints, double eps, std::vector& elts, std::vector& eltsIndex) const; - static void fillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt) throw(INTERP_KERNEL::Exception); - static void appendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector& ret); +/// @cond INTERNAL + typedef int (*DimM1DescNbrer)(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2); + MEDCouplingUMesh *buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const throw(INTERP_KERNEL::Exception); + static void FillInCompact3DMode(int spaceDim, int nbOfNodesInCell, const int *conn, const double *coo, double *zipFrmt) throw(INTERP_KERNEL::Exception); + static void AppendExtrudedCell(const int *connBg, const int *connEnd, int nbOfNodesPerLev, bool isQuad, std::vector& ret); + static void IntersectDescending2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, + std::vector< std::vector >& intersectEdge1, std::vector< std::vector >& subDiv2, + MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1, + 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 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); +/// @endcond private: //! this iterator stores current position in _nodal_connec array. mutable int _iterator; diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index a5da7e980..8e3ea8ff1 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -240,6 +240,7 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingUMesh::cellsByType; %newobject ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer; %newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity; +%newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity2; %newobject ParaMEDMEM::MEDCouplingUMesh::buildExtrudedMesh; %newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes; %newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshesOnSameCoords; @@ -1001,6 +1002,7 @@ namespace ParaMEDMEM void computeTypes() throw(INTERP_KERNEL::Exception); std::string reprConnectivityOfThis() const throw(INTERP_KERNEL::Exception); //tools + std::vector getQuadraticStatus() const throw(INTERP_KERNEL::Exception); DataArrayInt *findCellsIdsOnBoundary() const throw(INTERP_KERNEL::Exception); bool checkConsecutiveCellTypes() const throw(INTERP_KERNEL::Exception); DataArrayInt *rearrange2ConsecutiveCellTypes() throw(INTERP_KERNEL::Exception); @@ -1008,6 +1010,7 @@ namespace ParaMEDMEM DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception); DataArrayInt *zipConnectivityTraducer(int compType) throw(INTERP_KERNEL::Exception); MEDCouplingUMesh *buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception); + MEDCouplingUMesh *buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception); void orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception); bool isPresenceOfQuadratic() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *buildDirectionVectorField() const throw(INTERP_KERNEL::Exception); @@ -1259,6 +1262,26 @@ namespace ParaMEDMEM d3->incrRef(); return ret; } + + PyObject *buildDescendingConnectivity2() const throw(INTERP_KERNEL::Exception) + { + MEDCouplingAutoRefCountObjectPtr d0=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr d1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr d2=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr d3=DataArrayInt::New(); + MEDCouplingUMesh *m=self->buildDescendingConnectivity2(d0,d1,d2,d3); + PyObject *ret=PyTuple_New(5); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(m),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(d0),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(d1),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,3,SWIG_NewPointerObj(SWIG_as_voidptr(d2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,4,SWIG_NewPointerObj(SWIG_as_voidptr(d3),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + d0->incrRef(); + d1->incrRef(); + d2->incrRef(); + d3->incrRef(); + return ret; + } PyObject *emulateMEDMEMBDC(const MEDCouplingUMesh *nM1LevMesh) { @@ -1358,6 +1381,17 @@ namespace ParaMEDMEM return self->getCellIdsLyingOnNodes(da2->getConstPointer(),da2->getConstPointer()+da2->getNbOfElems(),fullyIn); } } + + static PyObject *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps) throw(INTERP_KERNEL::Exception) + { + DataArrayInt *cellNb1=0,*cellNb2=0; + MEDCouplingUMesh *mret=MEDCouplingUMesh::Intersect2DMeshes(m1,m2,eps,cellNb1,cellNb2); + PyObject *ret=PyTuple_New(3); + PyTuple_SetItem(ret,0,SWIG_NewPointerObj(SWIG_as_voidptr(mret),SWIGTYPE_p_ParaMEDMEM__MEDCouplingUMesh, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,1,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb1),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + PyTuple_SetItem(ret,2,SWIG_NewPointerObj(SWIG_as_voidptr(cellNb2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 )); + return ret; + } } void convertToPolyTypes(const std::vector& cellIdsToConvert) throw(INTERP_KERNEL::Exception); void convertAllToPoly(); -- 2.39.2