First working version. Still some doc and code re-factoring to be done.
}
/**
- * Reset the status of all edges (OUT, IN, ON) because they were potentially assignated
+ * Reset the status of all edges (OUT, IN, ON) because they were potentially assigned
* by the previous candidate processing.
*/
void ComposedEdge::initLocationsWithOther(const ComposedEdge& other) const
class ElementaryEdge;
class IteratorOnComposedEdge;
+ /**
+ * A set of quadratic or linear edges, described mainly by their connectivity
+ * The set is assumed to be connected, but not necessarily closed (i.e. not necessarily forming a closed polygon).
+ * Some methods however requires a closed form.
+ */
class ComposedEdge
{
friend class IteratorOnComposedEdge;
}
/*!
- * 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.
+ * This method builds from a start node, an end node and a direction a new ElementaryEdge.
*/
-ElementaryEdge *ElementaryEdge::BuildEdgeFromCrudeDataArray(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end)
+ElementaryEdge *ElementaryEdge::BuildEdgeFromStartEndDir(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end)
{
Edge *ptr=Edge::BuildEdgeFrom(start,end);
return new ElementaryEdge(ptr,direction);
std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const;
INTERPKERNEL_EXPORT void fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const;
- INTERPKERNEL_EXPORT static ElementaryEdge *BuildEdgeFromCrudeDataArray(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end);
+ INTERPKERNEL_EXPORT static ElementaryEdge *BuildEdgeFromStartEndDir(bool direction, INTERP_KERNEL::Node *start, INTERP_KERNEL::Node *end);
private:
bool _direction;
Edge *_ptr;
}
}
+void QuadraticPolygon::appendEdgeFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad,
+ const int *nodalBg, const double *coords,
+ const int segId, const std::vector<std::vector<int> >& intersectEdges)
+{
+ if(!isQuad)
+ {
+ bool direct=true;
+ int edgeId=segId;
+ const std::vector<int>& subEdge=intersectEdges[edgeId];
+ std::size_t nbOfSubEdges=subEdge.size()/2;
+ for(std::size_t j=0;j<nbOfSubEdges;j++)
+ appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
+ }
+ else
+ {
+ std::size_t nbOfSeg; //=std::distance(descBg,descEnd);
+ nbOfSeg = 1;
+ const double *st=coords+2*(nodalBg[0]);
+ INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
+ const double *endd=coords+2*(nodalBg[1]);
+ INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
+ const double *middle=coords+2*(nodalBg[2]);
+ INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
+ EdgeLin *e1,*e2;
+ e1=new EdgeLin(st0,middle0);
+ e2=new EdgeLin(middle0,endd0);
+ SegSegIntersector inters(*e1,*e2);
+ bool colinearity=inters.areColinears();
+ delete e1; delete e2;
+ //
+ bool direct=true;
+ int edgeId=segId;
+ const std::vector<int>& subEdge=intersectEdges[edgeId];
+ std::size_t nbOfSubEdges=subEdge.size()/2;
+ if(colinearity)
+ {
+ for(std::size_t j=0;j<nbOfSubEdges;j++)
+ appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
+ }
+ else
+ {
+ Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
+ for(std::size_t j=0;j<nbOfSubEdges;j++)
+ appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
+ e->decrRef();
+ }
+ st0->decrRef(); endd0->decrRef(); middle0->decrRef();
+ }
+}
+
+
void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector<int>& subEdge, const std::map<int,INTERP_KERNEL::Node *>& mapp)
{
std::size_t nbOfSubEdges=subEdge.size()/2;
{//it is not a quadratic subedge
Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
- ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
+ ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
pushBack(e);
}
else
}
/*!
- * 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
+ * This method builds a QP from the descending conn of a cell stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order, to give the
* orientation of edge.
+ * The information stored in intersectEdges1 and colinear1 is also used to know which intermediate nodes (intersection points with the other polygon) should be considered
+ * on the edges.
*/
void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
//appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
Node *start=(*mapp.find(idBg)).second;
Node *end=(*mapp.find(idEnd)).second;
- ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
+ ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
pushBack(e);
alreadyExistingIn2[descBg[i]].push_back(e);
}
}
}
+void QuadraticPolygon::buildFromCrudeDataArray3(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int segId, const std::vector<std::vector<int> >& intersectEdges,
+ const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
+ const std::vector< std::vector<int> >& colinear1,
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
+{
+ bool direct=true;
+ int edgeId=segId;//current edge id of pol2
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(segId);
+ if(it1!=alreadyExistingIn2.end())
+ {
+ const std::vector<INTERP_KERNEL::ElementaryEdge *>& edgesAlreadyBuilt=(*it1).second;
+ for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_iterator it3=edgesAlreadyBuilt.begin();it3!=edgesAlreadyBuilt.end();it3++)
+ {
+ Edge *ee=(*it3)->getPtr(); ee->incrRef();
+ pushBack(new ElementaryEdge(ee,(*it3)->getDirection()));
+ }
+ return;
+ }
+ bool directos=colinear1[edgeId].empty();
+ std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
+ int offset1=0;
+ if(!directos)
+ {// if the current edge of pol2 has one or more colinear edges part into pol1
+ const std::vector<int>& c=colinear1[edgeId];
+ std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
+ for(std::size_t j=0;j<nbOfEdgesIn1;j++)
+ {
+ int edgeId1=abs(descBg1[j])-1;
+ if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
+ {
+ idIns1.push_back(std::pair<int,std::pair<bool,int> >(edgeId1,std::pair<bool,int>(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2
+ //std::pair<edgeId1); direct1=descBg1[j]>0;
+ }
+ offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
+ }
+ directos=idIns1.empty();
+ }
+ if(directos)
+ {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
+ std::size_t oldSz=_sub_edges.size();
+ appendEdgeFromCrudeDataArray2(mapp,isQuad,nodalBg,coords,segId,intersectEdges);
+ std::size_t newSz=_sub_edges.size();
+ std::size_t zeSz=newSz-oldSz;
+ alreadyExistingIn2[segId].resize(zeSz);
+ std::list<ElementaryEdge *>::const_reverse_iterator it5=_sub_edges.rbegin();
+ for(std::size_t p=0;p<zeSz;p++,it5++)
+ alreadyExistingIn2[segId][zeSz-p-1]=*it5;
+ }
+ else
+ {//there is subpart of edge 'edgeId' of pol2 inside pol1
+ const std::vector<int>& subEdge=intersectEdges[edgeId];
+ std::size_t nbOfSubEdges=subEdge.size()/2;
+ for(std::size_t j=0;j<nbOfSubEdges;j++)
+ {
+ int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
+ int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
+ bool direction11,found=false;
+ bool direct1;//store if needed the direction in 1
+ int offset2;
+ std::size_t nbOfSubEdges1;
+ for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
+ {
+ int idIn1=(*it).first;//store if needed the cell id in 1
+ direct1=(*it).second.first;
+ offset1=(*it).second.second;
+ const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
+ nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
+ offset2=0;
+ for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
+ {//perform a loop on all subedges of pol1 that includes edge 'edgeId' of pol2. For the moment we iterate only on subedges of ['idIn1']... To improve
+ if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
+ { direction11=true; found=true; }
+ else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
+ { direction11=false; found=true; }
+ else
+ offset2++;
+ }
+ }
+ if(!found)
+ {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
+ //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
+ Node *start=(*mapp.find(idBg)).second;
+ Node *end=(*mapp.find(idEnd)).second;
+ ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
+ pushBack(e);
+ alreadyExistingIn2[segId].push_back(e);
+ }
+ else
+ {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
+ ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
+ Edge *ee=e->getPtr();
+ ee->incrRef();
+ ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
+ pushBack(e2);
+ alreadyExistingIn2[segId].push_back(e2);
+ }
+ }
+ }
+}
+
+/**
+ * @param edgeId current edge id of pol2
+ */
+void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray(const int edgeId, const std::vector<std::vector<int> >& intersectEdges,
+ const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
+ const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const
+{
+ std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
+ const std::vector<int>& c=colinear1[edgeId];
+ if(c.empty())
+ return;
+ const std::vector<int>& subEdge=intersectEdges[edgeId];
+ std::size_t nbOfSubEdges=subEdge.size()/2;
+ //
+ int offset1=0;
+ for(std::size_t j=0;j<nbOfEdgesIn1;j++)
+ {
+ int edgeId1=abs(descBg1[j])-1;
+ if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
+ {
+ for(std::size_t k=0;k<nbOfSubEdges;k++)
+ {
+ int idBg = subEdge[2*k];
+ int idEnd = subEdge[2*k+1];
+ int idIn1=edgeId1;
+ bool direct1=descBg1[j]>0;
+ const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
+ std::size_t nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
+ int offset2=0;
+ bool found=false;
+ for(std::size_t kk=0;kk<nbOfSubEdges1 && !found;kk++)
+ {
+ found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || \
+ (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
+ if(!found)
+ offset2++;
+ }
+ if(found)
+ {
+ ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
+ e->getPtr()->declareOn();
+ }
+ }
+ }
+ offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
+ }
+}
+
+
/*!
* Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
- * Method to find edges that are ON.
+ * Method to find edges that are ON, and mark them correspondingly on pol1.
*/
void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const
{
+ std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
std::size_t nbOfSeg=std::distance(descBg,descEnd);
for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
{
const std::vector<int>& subEdge=intersectEdges[edgeId];
std::size_t nbOfSubEdges=subEdge.size()/2;
//
- std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
int offset1=0;
for(std::size_t j=0;j<nbOfEdgesIn1;j++)
{
bool found=false;
for(std::size_t kk=0;kk<nbOfSubEdges1 && !found;kk++)
{
- found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
+ found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || \
+ (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
if(!found)
offset2++;
}
}
}
-void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, double xBary, double yBary, double fact, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const
+void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const
{
int nbOfNodesInPg=0;
bool presenceOfQuadratic=presenceOfQuadraticEdge();
for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
{
INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
- node->unApplySimilarity(xBary,yBary,fact);
+ //node->unApplySimilarity(xBary,yBary,fact);
addCoordsQuadratic.push_back((*node)[0]);
addCoordsQuadratic.push_back((*node)[1]);
conn.push_back(off+j);
}
/*!
- * 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 make the hypothesis that 'this' and 'other' are split 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'.
- * @param [in,out] edgesThis, parameter that keep informed the caller abount the edges in this not shared by the result of intersection of \a this with \a other
- * @param [in,out] edgesBoundaryOther, parameter that strores all edges in result of intersection that are not
+ * @param [in,out] edgesThis, parameter that keep informed the caller about the edges in this not shared by the result of intersection of \a this with \a other
+ * @param [in,out] edgesBoundaryOther, parameter that stores all edges in result of intersection that are not
*/
-void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther, const std::map<INTERP_KERNEL::Node *,int>& mapp, int idThis, int idOther, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nbThis, std::vector<int>& nbOther)
+void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis,
+ std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther, const std::map<INTERP_KERNEL::Node *,int>& mapp,
+ int idThis, int idOther, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn,
+ std::vector<int>& connI, std::vector<int>& nbThis, std::vector<int>& nbOther)
{
double xBaryBB, yBaryBB;
double fact=normalizeExt(&other, xBaryBB, yBaryBB);
- //Locate 'this' relative to 'other'
+ //Locate 'this' relative to 'other' (edges of 'this', aka 'pol1' are marked as IN or OUT)
other.performLocatingOperationSlow(*this); // without any assumption
- std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
+ std::vector<QuadraticPolygon *> res = BuildIntersectionPolygons(other,*this);
+ unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
+
for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
{
- (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
+ (*it)->appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI);
INTERP_KERNEL::IteratorOnComposedEdge it1(*it);
for(it1.first();!it1.finished();it1.next())
{
nbOther.push_back(idOther);
delete *it;
}
- unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
}
+void QuadraticPolygon::BuildPartitionFromZipList(const QuadraticPolygon & pol1, const std::list<INTERP_KERNEL::QuadraticPolygon *> & zip_list,
+ std::set<INTERP_KERNEL::Edge *> & edgesThis, std::set<INTERP_KERNEL::Edge *> & edgesBoundaryOther,
+ const std::map<INTERP_KERNEL::Node *,int>& mapp,
+ const int idThis, const std::vector<std::vector<int> > & mapZipTo2,
+ int offset, std::vector<double>& addCoordsQuadratic,
+ std::vector<int>& conn,std::vector<int>& connI,
+ std::vector<int>& nbThis, std::vector<int>& nbOther, std::vector<int>& nbOtherI)
+{
+ std::vector<QuadraticPolygon *> res;
+ if(!zip_list.empty())
+ ClosePolygonsSimple(pol1, zip_list, res, mapZipTo2, nbOther,nbOtherI);
+
+ for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
+ {
+ (*it)->appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI);
+ INTERP_KERNEL::IteratorOnComposedEdge it1(*it);
+ for(it1.first();!it1.finished();it1.next())
+ {
+ Edge *e=it1.current()->getPtr();
+ if(edgesThis.find(e)!=edgesThis.end())
+ edgesThis.erase(e);
+ else
+ {
+ if(edgesBoundaryOther.find(e)!=edgesBoundaryOther.end())
+ edgesBoundaryOther.erase(e);
+ else
+ edgesBoundaryOther.insert(e);
+ }
+ }
+ // Mapping to pol1 (mapping to pol2s has been handled in ClosePolygonsSimple).
+ nbThis.push_back(idThis);
+ delete (*it);
+ }
+}
+
+/** @sa ComposedEdge::initLocationsWithOther()
+ */
+void QuadraticPolygon::initLocationsWithSeveralOthers(const std::vector<QuadraticPolygon> & others) const
+{
+ std::set<Edge *> s1,s2;
+ for(std::list<ElementaryEdge *>::const_iterator it1=_sub_edges.begin();it1!=_sub_edges.end();it1++)
+ s1.insert((*it1)->getPtr());
+ std::vector<QuadraticPolygon>::const_iterator it_qp;
+ for(it_qp=others.begin(); it_qp!=others.end(); it_qp++)
+ for(std::list<ElementaryEdge *>::const_iterator it2=(*it_qp)._sub_edges.begin();it2!=(*it_qp)._sub_edges.end();it2++)
+ s2.insert((*it2)->getPtr());
+ initLocations();
+ for(it_qp=others.begin(); it_qp!=others.end(); it_qp++)
+ (*it_qp).initLocations();
+ std::vector<Edge *> s3;
+ std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::back_insert_iterator< std::vector<Edge *> >(s3));
+ for(std::vector<Edge *>::const_iterator it3=s3.begin();it3!=s3.end();it3++)
+ (*it3)->declareOn();
+}
+
+
/*!
* 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.
SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
//At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
performLocatingOperation(cpyOfOther);
- return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther);
+ return BuildIntersectionPolygons(cpyOfThis,cpyOfOther);
}
/*!
}
}
-void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const
+void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol1) const
{
- IteratorOnComposedEdge it(&pol2);
+ IteratorOnComposedEdge it(&pol1);
for(it.first();!it.finished();it.next())
{
ElementaryEdge *cur=it.current();
- cur->locateFullyMySelfAbsolute(*this);
+ cur->locateFullyMySelfAbsolute(*this); // *this=pol2=other
}
}
/*!
* Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned.
*
- * this : pol2 simplified.
* @param pol1 pol1 split.
* @param pol2 pol2 split.
*/
-std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
+std::vector<QuadraticPolygon *> QuadraticPolygon::BuildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2)
{
std::vector<QuadraticPolygon *> ret;
std::list<QuadraticPolygon *> pol2Zip=pol2.zipConsecutiveInSegments();
if(!pol2Zip.empty())
- closePolygons(pol2Zip,pol1,ret);
+ ClosePolygons(pol2Zip,pol1,pol2,ret);
else
{//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1
//do not overlap or pol1 is fully inside pol2. So in the first case no intersection, in the other case
//the intersection is pol1.
ElementaryEdge *e1FromPol1=pol1[0];
TypeOfEdgeLocInPolygon loc=FULL_ON_1;
- loc=e1FromPol1->locateFullyMySelf(*this,loc);
+ loc=e1FromPol1->locateFullyMySelf(pol2,loc);
if(loc==FULL_IN_1)
ret.push_back(new QuadraticPolygon(pol1));
}
return ret;
}
+std::list<QuadraticPolygon *> QuadraticPolygon::ZipConsecutiveSegments2(const std::vector<int> & candidates2,
+ const std::vector<QuadraticPolygon> & pol2s,
+ std::vector<std::vector<int> > & mapZipTo2)
+{
+ int prevI = -2;
+ TypeOfEdgeLocInPolygon prevLoc = FULL_UNKNOWN;
+ std::list<QuadraticPolygon *> ret;
+ std::vector<int>::const_iterator it2; int ii=0;
+ QuadraticPolygon *tmpQp = 0;
+ for(it2 = candidates2.begin(), ii = 0; it2 != candidates2.end(); it2++, ii++)
+ {
+ IteratorOnComposedEdge it_ii(const_cast<ComposedEdge*>(static_cast<const ComposedEdge*>(&pol2s[ii]))); // hmmm ... :-)
+ bool firstEdge = true;
+ for (; !it_ii.finished(); it_ii.next())
+ {
+ if(it_ii.current()->getLoc() == FULL_IN_1)
+ {
+ // only keep consecutive items from mesh2 (consecutive seg IDs and consecutive IN pieces)
+ if(prevLoc != FULL_IN_1 || (firstEdge && *it2 != prevI+1))
+ {
+ // the start node of a new QP has to be ON (otherwise we hit a piece of line starting in the middle of a cell)
+ if (it_ii.current()->getStartNode()->getLoc() == ON_1)
+ {
+ tmpQp = new QuadraticPolygon;
+ ret.push_back(tmpQp);
+ mapZipTo2.push_back(std::vector<int>(0));
+ }
+ else
+ continue; // firstEdge remains to its current value
+ }
+ tmpQp->pushBack(it_ii.current()->clone());
+ // Mapping - one quadratic polygon in ret is linked to several 1D cells:
+ if (std::find(mapZipTo2.back().begin(),mapZipTo2.back().end(), *it2) == mapZipTo2.back().end())
+ mapZipTo2.back().push_back(*it2);
+ }
+ prevLoc = it_ii.current()->getLoc();
+ firstEdge = false;
+ }
+ prevI = *it2;
+ }
+ // Now clean line pieces not ending ON:
+ for(std::list<QuadraticPolygon *>::iterator it3 = ret.begin(); it3 != ret.end(); it3++)
+ if ((*it3)->back()->getEndNode()->getLoc() != ON_1)
+ ret.erase(it3);
+
+ return ret;
+}
+
/*!
- * 'this' should be considered as pol2Simplified.
- * @param pol2zip is a list of set of edges (openned polygon) coming from split polygon 2.
+ * @param pol2zip is a list of set of edges (=an opened polygon) coming from split polygon 2.
* @param pol1 is split pol1.
* @param results the resulting \b CLOSED polygons.
*/
-void QuadraticPolygon::closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1,
- std::vector<QuadraticPolygon *>& results) const
+void QuadraticPolygon::ClosePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2,
+ std::vector<QuadraticPolygon *>& results)
{
bool directionKnownInPol1=false;
bool directionInPol1;
}
if(!directionKnownInPol1)
{
- if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1))
+ if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol1))
{ delete *iter; iter=pol2Zip.erase(iter); continue; }
else
directionKnownInPol1=true;
}
}
+void QuadraticPolygon::ClosePolygonsSimple(const QuadraticPolygon& pol1, const std::list<QuadraticPolygon *> & zip_list,
+ std::vector<QuadraticPolygon *>& results,const std::vector<std::vector<int> > & mapZipTo2,
+ std::vector<int>& nbOther, std::vector<int>& nbOtherI)
+{
+ // Register all edges of pol1 as 'alive' to start with. They will be consumed as they are
+ // added to the newly formed cells. Edges of pol2s (in zip_list) might be used more than once.
+ std::set<ElementaryEdge *> edges1_alive;
+ IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1));
+ for(; !it.finished(); it.next())
+ edges1_alive.insert(it.current());
+
+ int inf_loop_detect, i;
+ std::list<QuadraticPolygon *>::const_iterator itt;
+ for (inf_loop_detect = pol1.recursiveSize() + 1, itt = zip_list.begin(); itt != zip_list.end(); itt++)
+ inf_loop_detect += 2*(*itt)->recursiveSize();
+
+ for (i = 0, it.first(); !edges1_alive.empty() && i <= inf_loop_detect;)
+ {
+ QuadraticPolygon * qp = new QuadraticPolygon();
+ results.push_back(qp);
+ // Mapping initialization
+ int nbElemFrom2 = 0;
+ // Start anywhere valid (=alive) on pol1
+ ElementaryEdge * startE;
+ for (startE = it.current(); edges1_alive.find(it.current()) == edges1_alive.end(); it.nextLoop(), startE = it.current());
+
+ // Close the new cell:
+ do
+ {
+ ElementaryEdge *tmp_e = it.current();
+ qp->pushBack(tmp_e->clone());
+ edges1_alive.erase(tmp_e);
+ Node * nodeToTest = tmp_e->getEndNode();
+ it.nextLoop(); i++;
+ // Try to match the end node of the latest added edge from pol1 with one of the node from
+ // one of the ComposedEdge in zip_list
+ std::list<QuadraticPolygon *>::const_iterator ret;
+ bool direction;
+ int idx; // for mapping
+ ret = CheckInList2(nodeToTest, zip_list, direction, idx);
+ if (ret == zip_list.end())
+ continue;
+ size_t oldSz = nbOther.size();
+ size_t addSz = std::distance(mapZipTo2[idx].begin(), mapZipTo2[idx].end());
+ nbOther.resize(oldSz+addSz);
+ std::copy(mapZipTo2[idx].begin(), mapZipTo2[idx].end(), &nbOther[oldSz]);
+ nbElemFrom2 += addSz;
+ if (!direction)
+ (*ret)->reverse(); // TODO: discuss with Anthony - API of IteratorOnComposedEdge should allow an easy reverse looping?
+ // Switch to zip_list
+ IteratorOnComposedEdge it_pol2(*ret);
+ ElementaryEdge *tmp_e2;
+ for(/* it_pol2.first() */; !it_pol2.finished(); it_pol2.next(), i++)
+ {
+ tmp_e2 = it_pol2.current();
+ qp->pushBack(tmp_e2->clone());
+ }
+ // search where we landed on pol1 to resume looping on it.
+ if (!edges1_alive.empty())
+ {
+ int j;
+ for(j = 0; it.current()->getStartNode() != tmp_e2->getEndNode() && j < inf_loop_detect; j++) // better way of doing this?
+ it.nextLoop();
+ if (j >= inf_loop_detect) // could be a bit more precise and stop before ...
+ throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error - 1D intersection was looping infinitely!");
+ // sanity check - if we are not done with the cell completion, we should land
+ // on a live edge.
+ if(edges1_alive.find(it.current()) == edges1_alive.end() && it.current() != startE)
+ throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error ! Should never happen.");
+ }
+ }
+ while(it.current() != startE && !edges1_alive.empty() && i < inf_loop_detect);
+ // Complete mapping
+ nbOtherI.push_back(nbElemFrom2+nbOtherI.back());
+ }
+ if (i >= inf_loop_detect)
+ throw Exception("QuadraticPolygon::ClosePolygonsSimple() - Internal error - 1D intersection was looping infinitely!");
+}
+
/*!
* 'this' is expected to be set of edges (not closed) of pol2 split.
*/
-bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
+bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
{
IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
bool found=false;
it.next();
}
if(!found)
- throw Exception("Internal error : polygons uncompatible each others. Should never happend");
+ throw Exception("Internal error: polygons incompatible with each others. Should never happen!");
//Ok we found correspondance between this and pol1. Searching for right direction to close polygon.
ElementaryEdge *e=_sub_edges.back();
if(e->getLoc()==FULL_ON_1)
}
/*!
- * This method fills as much as possible 'this' (part of pol2 split) with edges of 'pol1Splitted'.
+ * This method fills as much as possible 'this' (a sub-part of pol2 split) with edges of 'pol1Splitted'.
*/
std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
std::list<QuadraticPolygon *>::iterator iStart,
return iEnd;
}
+std::list<QuadraticPolygon *>::const_iterator QuadraticPolygon::CheckInList2(Node *n, const std::list<QuadraticPolygon *> & zip_list,
+ bool & direction, int & index)
+{
+ direction = true;
+ index = 0;
+ for(std::list<QuadraticPolygon *>::const_iterator iter = zip_list.begin(); iter!=zip_list.end(); iter++, index++)
+ {
+ if ((*iter)->front()->getStartNode() == n)
+ return iter;
+ else
+ {
+ if ((*iter)->back()->getEndNode() == n)
+ {
+ direction = false;
+ return iter;
+ }
+ }
+ }
+ return zip_list.end();
+}
+
+
void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary, const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, int idThis,
std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2)
{
std::list<QuadraticPolygon *> pol1Zip;
if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty())
{
- pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
+ pol1.appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
return ;
}
while(!notUsedInPol1L.empty())
{
if((*it1)->getStartNode()==(*it1)->getEndNode())
{
- (*it1)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
+ (*it1)->appendCrudeData(mapp,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
delete *it6;
delete *it1;
const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
const std::vector< std::vector<int> >& colinear1,
std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2);
+ INTERPKERNEL_EXPORT void buildFromCrudeDataArray3(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int segId, const std::vector<std::vector<int> >& intersectEdges,
+ const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
+ const std::vector< std::vector<int> >& colinear1,
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2);
+ INTERPKERNEL_EXPORT void updateLocOfEdgeFromCrudeDataArray(const int edgeId, const std::vector<std::vector<int> >& intersectEdges,
+ const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
+ const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const;
INTERPKERNEL_EXPORT void updateLocOfEdgeFromCrudeDataArray2(const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges, const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const;
INTERPKERNEL_EXPORT void appendEdgeFromCrudeDataArray(std::size_t edgeId, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges);
+ INTERPKERNEL_EXPORT void appendEdgeFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad,
+ const int *nodalBg, const double *coords,
+ const int segId, const std::vector<std::vector<int> >& intersectEdges);
INTERPKERNEL_EXPORT void appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector<int>& subEdge, const std::map<int,INTERP_KERNEL::Node *>& mapp);
- INTERPKERNEL_EXPORT void appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, double xBary, double yBary, double fact, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const;
+ INTERPKERNEL_EXPORT void appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const;
INTERPKERNEL_EXPORT void buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther, const std::map<INTERP_KERNEL::Node *,int>& mapp, int idThis, int idOther, int offset,
std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2);
+ INTERPKERNEL_EXPORT static void BuildPartitionFromZipList(const QuadraticPolygon & pol1, const std::list<INTERP_KERNEL::QuadraticPolygon *> & zip_list,
+ std::set<INTERP_KERNEL::Edge *> & edgesThis, std::set<INTERP_KERNEL::Edge *> & edgesBoundaryOther,
+ const std::map<INTERP_KERNEL::Node *,int>& mapp,
+ const int idThis, const std::vector<std::vector<int> > & mapZipTo2,
+ int offset, std::vector<double>& addCoordsQuadratic,
+ std::vector<int>& conn,std::vector<int>& connI,
+ std::vector<int>& nbThis, std::vector<int>& nbOther, std::vector<int>& nbOtherI);
+ INTERPKERNEL_EXPORT void initLocationsWithSeveralOthers(const std::vector<QuadraticPolygon> & others) const;
//
INTERPKERNEL_EXPORT double intersectWith(const QuadraticPolygon& other) const;
INTERPKERNEL_EXPORT double intersectWith(const QuadraticPolygon& other, double* barycenter) const;
INTERPKERNEL_EXPORT void performLocatingOperation(QuadraticPolygon& pol2) const;
INTERPKERNEL_EXPORT void performLocatingOperationSlow(QuadraticPolygon& pol2) const;
INTERPKERNEL_EXPORT static void SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits);
- INTERPKERNEL_EXPORT std::vector<QuadraticPolygon *> buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const;
- INTERPKERNEL_EXPORT bool amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction);
+ INTERPKERNEL_EXPORT static std::vector<QuadraticPolygon *> BuildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2);
+ INTERPKERNEL_EXPORT bool haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction);
INTERPKERNEL_EXPORT static void ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary, const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, int idThis,
std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2);
+ INTERPKERNEL_EXPORT static std::list<QuadraticPolygon *> ZipConsecutiveSegments2(const std::vector<int> & candidates2, const std::vector<QuadraticPolygon> & pol2s,
+ std::vector<std::vector<int> > & mapZipTo2);
protected:
std::list<QuadraticPolygon *> zipConsecutiveInSegments() const;
void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const;
- void closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1, std::vector<QuadraticPolygon *>& results) const;
+ static void ClosePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2, std::vector<QuadraticPolygon *>& results);
+ static void ClosePolygonsSimple(const QuadraticPolygon& pol1, const std::list<INTERP_KERNEL::QuadraticPolygon *> & zip_list,
+ std::vector<QuadraticPolygon *>& results, const std::vector<std::vector<int> > & mapZipTo2,
+ std::vector<int>& nbOther, std::vector<int>& nbOtherI);
template<class EDGES>
static void UpdateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2,
const EDGES *e1, const EDGES *e2);
+
std::list<QuadraticPolygon *>::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
std::list<QuadraticPolygon *>::iterator iStart,
std::list<QuadraticPolygon *>::iterator iEnd,
bool direction);
static std::list<QuadraticPolygon *>::iterator CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
std::list<QuadraticPolygon *>::iterator iEnd);
+ static std::list<QuadraticPolygon *>::const_iterator CheckInList2(Node *n, const std::list<QuadraticPolygon *> & zip_list,
+ bool & direction, int & index);
};
}
}
/**
- * Construct a mapping between set of Nodes and the standart MEDCoupling connectivity format (c, cI).
+ * Construct a mapping between set of Nodes and the standard MEDCoupling connectivity format (c, cI).
*/
void MEDCouplingUMeshBuildQPFromMesh3(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
const int *desc1Bg, const int *desc1End, const std::vector<std::vector<int> >& intesctEdges1,
}
}
}
+
+ /**
+ * Same as above for a 1D cell made of a single segment
+ * TODO: factorize with the above
+ */
+ void MEDCouplingUMeshBuildQPFromMesh4(const double *coo1, int offset1, const double *coo2, int offset2, const std::vector<double>& addCoo,
+ const int segId, const std::vector<std::vector<int> >& intesctEdges1,
+ /*output*/std::map<INTERP_KERNEL::Node *,int>& mapp, std::map<int,INTERP_KERNEL::Node *>& mappRev)
+ {
+ for(std::vector<int>::const_iterator it1=intesctEdges1[segId].begin();it1!=intesctEdges1[segId].end();it1++)
+ {
+ std::map<int,INTERP_KERNEL::Node *>::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;
+ }
+ }
+ }
+
}
/// @endcond
MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Intersect2D",2);
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
- MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
ret->setConnectivity(conn,connI,true);
ret->setCoords(coo);
return ret.retn();
}
+MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
+ double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2)
+{
+ m1->checkFullyDefined();
+ m2->checkFullyDefined();
+ if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=1 || m2->getSpaceDimension()!=2)
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine works on unstructured meshes m1 with (meshdim, spacedim) = (2,2) and m2 with (meshdim, spacedim) = (1,2)!");
+
+ // Check that m2 is a line (and not a more complex 1D mesh) - each point is used at most by 2 segments:
+ DataArrayInt * d = DataArrayInt::New(), * dI = DataArrayInt::New();
+ DataArrayInt *rD = DataArrayInt::New(), * rDI = DataArrayInt::New();
+ MEDCouplingUMesh * m_points;
+ m_points = m2->buildDescendingConnectivity(d, dI, rD, rDI);
+ DataArrayInt * dsi = rDI->deltaShiftIndex();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dsii = dsi->getIdsNotInRange(0,3);
+ m_points->decrRef(); d->decrRef(); dI->decrRef(); rD->decrRef(); rDI->decrRef(); dsi->decrRef();
+ if (dsii->getNumberOfTuples())
+ throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshWith1DLine only work with mesh2 being a (piecewise) connected line!");
+
+ // Step 1: compute all edge intersections (new nodes)
+ std::vector< std::vector<int> > intersectEdge1, colinear2, subDiv2;
+ MEDCouplingUMesh *m1Desc=0;
+ DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0;
+ std::vector<double> addCoo,addCoordsQuadratic; // coordinates of newly created nodes
+ INTERP_KERNEL::QUADRATIC_PLANAR::_precision=eps;
+ INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=eps;
+ IntersectDescending2DMeshWith1DMesh(m1,m2,eps,intersectEdge1,colinear2, subDiv2,
+ m1Desc,desc1,descIndx1,revDesc1,revDescIndx1,
+ addCoo);
+ revDesc1->decrRef(); revDescIndx1->decrRef();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd5(m1Desc);
+
+ // Step 2: re-order newly created nodes according to the ordering found in m2
+ std::vector< std::vector<int> > intersectEdge2;
+ BuildIntersectEdges(m1Desc,m2,addCoo,subDiv2,intersectEdge2);
+ subDiv2.clear(); dd5=0;
+
+ // Step 3: this is where we have a significant difference with Intersect2DMeshes
+ std::vector<int> cr,crI; //no DataArrayInt because interface with Geometric2D
+ std::vector<int> cNb1,cNb2,cNbI2; //no DataArrayInt because interface with Geometric2D
+ cNbI2.push_back(0);
+ Build2DCellsFrom1DCut(eps,m1,desc1->getConstPointer(),descIndx1->getConstPointer(),intersectEdge1,
+ colinear2,m2,intersectEdge2,addCoo,
+ /* outputs -> */addCoordsQuadratic,cr,crI,cNb1,cNb2,cNbI2);
+
+ // Step 4: Prepare final result:
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCooDa=DataArrayDouble::New();
+ addCooDa->alloc((int)(addCoo.size())/2,2);
+ std::copy(addCoo.begin(),addCoo.end(),addCooDa->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> addCoordsQuadraticDa=DataArrayDouble::New();
+ addCoordsQuadraticDa->alloc((int)(addCoordsQuadratic.size())/2,2);
+ std::copy(addCoordsQuadratic.begin(),addCoordsQuadratic.end(),addCoordsQuadraticDa->getPointer());
+ std::vector<const DataArrayDouble *> coordss(4);
+ coordss[0]=m1->getCoords(); coordss[1]=m2->getCoords(); coordss[2]=addCooDa; coordss[3]=addCoordsQuadraticDa;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::Aggregate(coordss);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> ret=MEDCouplingUMesh::New("Intersect2D",2);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> conn=DataArrayInt::New(); conn->alloc((int)cr.size(),1); std::copy(cr.begin(),cr.end(),conn->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> connI=DataArrayInt::New(); connI->alloc((int)crI.size(),1); std::copy(crI.begin(),crI.end(),connI->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c1=DataArrayInt::New(); c1->alloc((int)cNb1.size(),1); std::copy(cNb1.begin(),cNb1.end(),c1->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> c2=DataArrayInt::New(); c2->alloc((int)cNb2.size(),1); std::copy(cNb2.begin(),cNb2.end(),c2->getPointer());
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> cI2=DataArrayInt::New(); cI2->alloc((int)cNbI2.size(),1); std::copy(cNbI2.begin(),cNbI2.end(),cI2->getPointer());
+ ret->setConnectivity(conn,connI,true);
+ ret->setCoords(coo);
+ cellNb1=c1.retn(); cellNb2=c2.retn();cellNbI2=cI2.retn();
+ return ret.retn();
+}
/**
* Private. Third step of the partitioning algorithm (Intersect2DMeshes): reconstruct full 2D cells from the
const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
// Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,/* output */mapp,mappRev);
- // pol1 is the full cell from mesh2, in QP format, with all the additional intersecting nodes.
+ // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
//
INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
// Complete mapping with elements coming from the current cell it2 in mesh2:
- MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
+ MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc2+descIndx2[*it2],
+ desc2+descIndx2[*it2+1],intesctEdges2,/* output */mapp,mappRev);
// pol2 is the new QP in the final merged result.
pol2s[ii].buildFromCrudeDataArray2(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,desc2+descIndx2[*it2],desc2+descIndx2[*it2+1],intesctEdges2,
pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
}
}
+void MEDCouplingUMesh::Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
+ const MEDCouplingUMesh *m2, const std::vector<std::vector<int> >& intesctEdges2,
+ const std::vector<double>& addCoords,
+ std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI,
+ std::vector<int>& cNb1, std::vector<int>& cNb2, std::vector<int>& cNbI2)
+{
+ static const int SPACEDIM=2;
+ 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();
+ const int *conn2=m2->getNodalConnectivity()->getConstPointer();
+ const int *connI2=m2->getNodalConnectivityIndex()->getConstPointer();
+ int offset2=offset1+m2->getNumberOfNodes();
+ int offset3=offset2+((int)addCoords.size())/2;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
+ const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
+ // Here a BBTree on 1D-segments from the tool mesh, not on segments:
+ BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2->getNumberOfCells(),eps);
+ int ncell1=m1->getNumberOfCells();
+ crI.push_back(0);
+ // for all 2D cells in base mesh, find intersecting segments
+ for(int i=0;i<ncell1;i++)
+ {
+ std::vector<int> candidates2;
+ myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
+ // Sort the candidates: this will help reconstructing connected lines traversing cell(i) in m1.
+ std::sort( candidates2.begin(), candidates2.end() );
+ std::map<INTERP_KERNEL::Node *,int> mapp;
+ std::map<int,INTERP_KERNEL::Node *> mappRev;
+ // conversion to the geometric DS format:
+ INTERP_KERNEL::QuadraticPolygon pol1;
+ INTERP_KERNEL::NormalizedCellType typ=(INTERP_KERNEL::NormalizedCellType)conn1[connI1[i]];
+ const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(typ);
+ // Populate mapp and mappRev with nodes from the current cell (i) from mesh1 - this also builds the Node* objects:
+ MEDCouplingUMeshBuildQPFromMesh3(coo1,offset1,coo2,offset2,addCoords,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,
+ /* output-->*/mapp,mappRev);
+ // pol1 is the full cell from mesh1, in QP format, with all the additional intersecting nodes.
+ pol1.buildFromCrudeDataArray(mappRev,cm.isQuadratic(),conn1+connI1[i]+1,coo1,
+ desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1);
+ //
+ std::set<INTERP_KERNEL::Edge *> edges1;// store all edges of pol1 that are NOT consumed by intersect cells.
+ // If any after the for-loop over candidates2 -> it remains unhandled parts of pol1 which should appear in the result (sse computeResidual below)
+ std::set<INTERP_KERNEL::Edge *> edgesBoundary2;// store all edges that are on boundary of (pol2 intersect pol1) minus edges on pol1.
+ INTERP_KERNEL::IteratorOnComposedEdge it1(&pol1);
+ for(it1.first();!it1.finished();it1.next())
+ edges1.insert(it1.current()->getPtr());
+ //
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> > edgesIn2ForShare; // common edges
+ std::vector<INTERP_KERNEL::QuadraticPolygon> pol2s(candidates2.size());
+ std::vector<int>::iterator it2; int ii=0;
+ for(it2=candidates2.begin();it2!=candidates2.end();it2++,ii++)
+ {
+ INTERP_KERNEL::NormalizedCellType typ2=(INTERP_KERNEL::NormalizedCellType)conn2[connI2[*it2]];
+ const INTERP_KERNEL::CellModel& cm2=INTERP_KERNEL::CellModel::GetCellModel(typ2);
+ // Complete mapping with elements coming from the current segment it2 in linear mesh2:
+ MEDCouplingUMeshBuildQPFromMesh4(coo1,offset1,coo2,offset2,addCoords,*it2,intesctEdges2,
+ /* output */mapp,mappRev);
+ // pol2 is the new QP in the final merged result.
+ pol2s[ii].buildFromCrudeDataArray3(mappRev,cm2.isQuadratic(),conn2+connI2[*it2]+1,coo2,*it2,intesctEdges2,
+ pol1,desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2, /* output */ edgesIn2ForShare);
+ }
+ pol1.initLocationsWithSeveralOthers(pol2s);
+ for(it2 = candidates2.begin(), ii = 0; it2 != candidates2.end(); it2++, ii++)
+ {
+ pol2s[ii].updateLocOfEdgeFromCrudeDataArray(*it2 ,intesctEdges2,pol1, desc1+descIndx1[i],desc1+descIndx1[i+1],intesctEdges1,colinear2);
+ double xBaryBB, yBaryBB;
+ double fact=pol1.normalizeExt(&pol2s[ii], xBaryBB, yBaryBB);
+ //Locate pol2s[ii] relative to pol1 (this is done the other way around in the 2D/2D intersection)
+ pol1.performLocatingOperationSlow(pol2s[ii]);
+ pol1.unApplyGlobalSimilarityExt(pol2s[ii],xBaryBB,yBaryBB,fact);
+ }
+ // At this point all edges in (all elements of) pol2s are fully localized.
+ // Zip consecutive IN segments from the list of candidates, and discard the pieces which are ending nowhere (e.g. an IN polyline
+ // with an end-point in the middle of the cell).
+ std::vector<std::vector<int> > mapZipTo2;
+ std::list<INTERP_KERNEL::QuadraticPolygon *> zip_list = INTERP_KERNEL::QuadraticPolygon::ZipConsecutiveSegments2(candidates2, pol2s, mapZipTo2);
+
+ // Now the core of the algo - main output is in cr, crI, cNb1, cNb2 and cNbI2.
+ INTERP_KERNEL::QuadraticPolygon::BuildPartitionFromZipList(pol1, zip_list, edges1, edgesBoundary2, mapp,
+ i, mapZipTo2,
+ offset3,addCoordsQuadratic,cr,crI, cNb1,cNb2, cNbI2);
+ // Deals with remaining (non-consumed) edges from m1: these are the edges that were never touched
+ // by m2 but that we still want to keep in the final result.
+ if(!edges1.empty())
+ {
+ size_t oldSz = cNb2.size();
+ try
+ {
+ INTERP_KERNEL::QuadraticPolygon::ComputeResidual(pol1,edges1,edgesBoundary2,mapp,offset3,i,addCoordsQuadratic,cr,crI,cNb1,cNb2);
+ }
+ catch(INTERP_KERNEL::Exception& e)
+ {
+ std::ostringstream oss; oss << "Error when computing residual of cell #" << i << " in source/m1 mesh ! Maybe the neighbours of this cell in mesh are not well connected !\n" << "The deep reason is the following : " << e.what();
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ // ComputeResidual has pushed some -1 at the back of cNb2. Replace this with void entries in cNbI2:
+ // TODO: discuss the format with Anthony.
+ size_t newSz = cNb2.size();
+ cNb2.resize(oldSz);
+ int val = (oldSz == 0) ? 0 : cNbI2.back();
+ for(size_t sss=0; sss<newSz-oldSz; sss++, cNbI2.push_back(val));
+ }
+ // Memory clean-up
+ for(std::map<int,INTERP_KERNEL::Node *>::const_iterator it=mappRev.begin();it!=mappRev.end();it++)
+ (*it).second->decrRef();
+ for (std::list<INTERP_KERNEL::QuadraticPolygon *>::const_iterator itt = zip_list.begin(); itt != zip_list.end(); itt++)
+ delete(*itt);
+ }
+}
+
+
void IKGeo2DInternalMapper2(INTERP_KERNEL::Node *n, const std::map<INTERP_KERNEL::Node *,int>& m, int forbVal0, int forbVal1, std::vector<int>& isect)
{
std::map<INTERP_KERNEL::Node *,int>::const_iterator it(m.find(n));
intersectEdge1.resize(nDescCell1);
colinear2.resize(nDescCell2);
subDiv2.resize(nDescCell2);
- BBTree<SPACEDIM,int> myTree(bbox2,0,0,m2Desc->getNumberOfCells(),-eps);
+ BBTree<SPACEDIM,int> myTree(bbox2,0,0,nDescCell2,-eps);
std::vector<int> candidates1(1);
int offset1=m1->getNumberOfNodes();
m2Desc->incrRef(); desc2->incrRef(); descIndx2->incrRef(); revDesc2->incrRef(); revDescIndx2->incrRef();
}
+void MEDCouplingUMesh::IntersectDescending2DMeshWith1DMesh(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
+ std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
+ MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
+ std::vector<double>& addCoo)
+{
+ static const int SPACEDIM=2;
+ // Build desc connectivity of first mesh only
+ desc1=DataArrayInt::New(); descIndx1=DataArrayInt::New();
+ revDesc1=DataArrayInt::New(); revDescIndx1=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> dd1(desc1),dd2(descIndx1),dd3(revDesc1),dd4(revDescIndx1);
+ m1Desc=m1->buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> dd9(m1Desc);
+ const int *c1=m1Desc->getNodalConnectivity()->getConstPointer();
+ const int *ci1=m1Desc->getNodalConnectivityIndex()->getConstPointer();
+
+ // Build BB tree of all edges in the tool mesh (second mesh)
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> bbox1Arr(m1Desc->getBoundingBoxForBBTree()),bbox2Arr(m2->getBoundingBoxForBBTree());
+ const double *bbox1(bbox1Arr->begin()),*bbox2(bbox2Arr->begin());
+ int nDescCell1=m1Desc->getNumberOfCells();
+ int nCell2=m2->getNumberOfCells();
+ intersectEdge1.resize(nDescCell1);
+ colinear2.resize(nCell2);
+ subDiv2.resize(nCell2);
+ BBTree<SPACEDIM,int> myTree(bbox2,0,0,nCell2,-eps);
+
+ std::vector<int> candidates1(1);
+ int offset1=m1->getNumberOfNodes();
+ int offset2=offset1+m2->getNumberOfNodes();
+ for(int i=0;i<nDescCell1;i++) // for all edges in the first mesh
+ {
+ std::vector<int> candidates2; // edges of mesh2 candidate for intersection
+ myTree.getIntersectingElems(bbox1+i*2*SPACEDIM,candidates2);
+ if(!candidates2.empty()) // candidates2 holds edges from the second mesh potentially intersecting current edge i in mesh1
+ {
+ std::map<INTERP_KERNEL::Node *,int> map1,map2;
+ // pol2 is not necessarily a closed polygon: just a set of (quadratic) edges (same as candidates2) in the Geometric DS format
+ INTERP_KERNEL::QuadraticPolygon *pol2=MEDCouplingUMeshBuildQPFromMesh(m2,candidates2,map2);
+ candidates1[0]=i;
+ INTERP_KERNEL::QuadraticPolygon *pol1=MEDCouplingUMeshBuildQPFromMesh(m1Desc,candidates1,map1);
+ // This following part is to avoid that some removed nodes (for example due to a merge between pol1 and pol2) are replaced by a newly created one
+ // This trick guarantees that Node * are discriminant (i.e. form a unique identifier)
+ std::set<INTERP_KERNEL::Node *> nodes;
+ pol1->getAllNodes(nodes); pol2->getAllNodes(nodes);
+ std::size_t szz(nodes.size());
+ std::vector< MEDCouplingAutoRefCountObjectPtr<INTERP_KERNEL::Node> > nodesSafe(szz);
+ std::set<INTERP_KERNEL::Node *>::const_iterator itt(nodes.begin());
+ for(std::size_t iii=0;iii<szz;iii++,itt++)
+ { (*itt)->incrRef(); nodesSafe[iii]=*itt; }
+ // end of protection
+ // Performs egde cutting:
+ pol1->splitAbs(*pol2,map1,map2,offset1,offset2,candidates2,intersectEdge1[i],i,colinear2,subDiv2,addCoo);
+ delete pol2;
+ delete pol1;
+ }
+ else
+ intersectEdge1[i].insert(intersectEdge1[i].end(),c1+ci1[i]+1,c1+ci1[i+1]);
+ }
+ m1Desc->incrRef(); desc1->incrRef(); descIndx1->incrRef(); revDesc1->incrRef(); revDescIndx1->incrRef();
+}
+
+
/*!
* This method performs the 2nd step of Partition of 2D mesh.
* This method has 4 inputs :
MEDCOUPLING_EXPORT static void ComputeVecAndPtOfFace(double eps, const double *coords, const int *begin, const int *end, double *v, double *p);
MEDCOUPLING_EXPORT static void TryToCorrectPolyhedronOrientation(int *begin, int *end, const double *coords);
MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2);
+ MEDCOUPLING_EXPORT static MEDCouplingUMesh *Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2,
+ double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2, DataArrayInt *&cellNbI2);
MEDCOUPLING_EXPORT static bool BuildConvexEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, DataArrayInt *nodalConnecOut);
MEDCOUPLING_EXPORT static bool RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, const int *idsToRemoveEnd, DataArrayInt *arr, DataArrayInt *arrIndx, int offsetForRemoval=0);
MEDCOUPLING_EXPORT static void ExtractFromIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
std::vector<double>& addCoo,
MEDCouplingUMesh *& m2Desc, DataArrayInt *&desc2, DataArrayInt *&descIndx2, DataArrayInt *&revDesc2, DataArrayInt *&revDescIndx2);
+ static void IntersectDescending2DMeshWith1DMesh(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
+ std::vector< std::vector<int> >& intersectEdge1, std::vector< std::vector<int> >& colinear2, std::vector< std::vector<int> >& subDiv2,
+ MEDCouplingUMesh *& m1Desc, DataArrayInt *&desc1, DataArrayInt *&descIndx1, DataArrayInt *&revDesc1, DataArrayInt *&revDescIndx1,
+ std::vector<double>& addCoo);
static void BuildIntersectEdges(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, const std::vector<double>& addCoo, const std::vector< std::vector<int> >& subDiv, std::vector< std::vector<int> >& intersectEdge);
static void BuildIntersecting2DCellsFromEdges(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
const MEDCouplingUMesh *m2, const int *desc2, const int *descIndx2, const std::vector<std::vector<int> >& intesctEdges2,
const std::vector<double>& addCoords,
std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI, std::vector<int>& cNb1, std::vector<int>& cNb2);
+ static void Build2DCellsFrom1DCut(double eps, const MEDCouplingUMesh *m1, const int *desc1, const int *descIndx1, const std::vector<std::vector<int> >& intesctEdges1, const std::vector< std::vector<int> >& colinear2,
+ const MEDCouplingUMesh *m2, const std::vector<std::vector<int> >& intesctEdges2,
+ const std::vector<double>& addCoords,
+ std::vector<double>& addCoordsQuadratic, std::vector<int>& cr, std::vector<int>& crI,
+ std::vector<int>& cNb1, std::vector<int>& cNb2, std::vector<int>& cNbI2);
static void AssemblyForSplitFrom3DCurve(const std::vector<int>& cut3DCurve, std::vector<int>& nodesOnPlane, const int *nodal3DSurf, const int *nodalIndx3DSurf,
const int *nodal3DCurve, const int *nodalIndx3DCurve,
const int *desc, const int *descIndx, std::vector< std::pair<int,int> >& cut3DSurf);
//
m->decrRef();
}
+
+/**
+ * Two square cells intersected by a vertical line crossing only one cell.
+ */
+void MEDCouplingBasicsTest5::testIntersect2DMeshWith1DLine1()
+{
+ MEDCouplingCMesh *m1c = MEDCouplingCMesh::New();
+ DataArrayDouble *coordX = DataArrayDouble::New();
+ const double arrX[3] = {-1., 1., 2};
+ coordX->alloc(3,1);
+ std::copy(arrX,arrX+3,coordX->getPointer());
+ m1c->setCoordsAt(0,coordX);
+ DataArrayDouble *coordY = DataArrayDouble::New();
+ const double arrY[4] = {0., 2.};
+ coordY->alloc(2,1);
+ std::copy(arrY,arrY+2,coordY->getPointer());
+ m1c->setCoordsAt(1,coordY);
+ MEDCouplingUMesh *m1 = m1c->buildUnstructured();
+
+ // A simple line:
+ MEDCouplingUMesh *m2 = MEDCouplingUMesh::New("bla", 1);
+ const int conn2A[7] = {INTERP_KERNEL::NORM_SEG2,0,1,INTERP_KERNEL::NORM_SEG3,1,2,3};
+ const int connI2A[3] = {0,3,7};
+ const double coo2A[8] = {0.,-1.0, 0.,1., 0.,3., 0.5,2.2};
+ DataArrayDouble *coord2 = DataArrayDouble::New();
+ DataArrayInt *conn2 = DataArrayInt::New();
+ DataArrayInt *connI2 = DataArrayInt::New();
+ coord2->alloc(4,2); conn2->alloc(7,1); connI2->alloc(3, 1);
+ std::copy(coo2A ,coo2A+8, coord2->getPointer());
+ std::copy(conn2A ,conn2A+7, conn2->getPointer());
+ std::copy(connI2A,connI2A+3, connI2->getPointer());
+ m2->setCoords(coord2);
+ m2->setConnectivity(conn2, connI2);
+
+ // End of construction of input meshes m1bis and m2 -> start of specific part of the test
+ DataArrayInt * map1 = 0, * map2 = 0, * mapI2 = 0;
+ MEDCouplingUMesh *m3 = MEDCouplingUMesh::Intersect2DMeshWith1DLine(m1, m2, 1e-10, map1, map2, mapI2);
+ bool areNodesMerged; int newNbNodes;
+ DataArrayInt * d_tmp = m3->mergeNodes(1.0e-8, areNodesMerged, newNbNodes);
+ d_tmp->decrRef();
+ CPPUNIT_ASSERT_EQUAL(3,m3->getNumberOfCells());
+ CPPUNIT_ASSERT_EQUAL(20,m3->getNumberOfNodes());
+ CPPUNIT_ASSERT_EQUAL(2,m3->getSpaceDimension());
+
+ const int exp1[3] = {0,0,1};
+ const int exp2[4] = {0,1,0,1};
+ const int expI2[4] = {0,2,4,4};
+ CPPUNIT_ASSERT_EQUAL(3, map1->getNumberOfTuples());
+ CPPUNIT_ASSERT_EQUAL(4, map2->getNumberOfTuples());
+ CPPUNIT_ASSERT_EQUAL(4, mapI2->getNumberOfTuples());
+ CPPUNIT_ASSERT(std::equal(exp1,exp1+3, map1->getConstPointer()));
+ CPPUNIT_ASSERT(std::equal(exp2,exp2+4, map2->getConstPointer()));
+ CPPUNIT_ASSERT(std::equal(expI2,expI2+4, mapI2->getConstPointer()));
+
+ const int expConn[27] = {32, 1, 10, 7, 11, 4, 12, 13, 14, 15, 16, 32, 10, 0, 3, 11, 7, 17, 18, 19, 14, 13, 5, 2, 1, 4, 5};
+ const int expConnI[4] = {0, 11, 22, 27};
+ CPPUNIT_ASSERT_EQUAL(27, m3->getNodalConnectivity()->getNumberOfTuples());
+ CPPUNIT_ASSERT_EQUAL(4, m3->getNodalConnectivityIndex()->getNumberOfTuples());
+ CPPUNIT_ASSERT(std::equal(expConn,expConn+27, m3->getNodalConnectivity()->getConstPointer()));
+ CPPUNIT_ASSERT(std::equal(expConnI,expConnI+4, m3->getNodalConnectivityIndex()->getConstPointer()));
+
+ map1->decrRef();
+ map2->decrRef();
+ mapI2->decrRef();
+ m3->decrRef();
+ //
+ m2->decrRef();
+ m1->decrRef();
+ coordX->decrRef();
+ coordY->decrRef();
+ m1c->decrRef();
+ coord2->decrRef();
+ conn2->decrRef();
+ connI2->decrRef();
+}
+
CPPUNIT_TEST( testDAIBuildSubstractionOptimized1 );
CPPUNIT_TEST( testDAIIsStrictlyMonotonic1 );
CPPUNIT_TEST( testSimplexize3 );
+ CPPUNIT_TEST( testIntersect2DMeshWith1DLine1 );
CPPUNIT_TEST_SUITE_END();
public:
void testUMeshTessellate2D1();
void testDAIBuildSubstractionOptimized1();
void testDAIIsStrictlyMonotonic1();
void testSimplexize3();
+ void testIntersect2DMeshWith1DLine1();
};
}
self.assertTrue(m.getNodalConnectivityIndex().isEqual(DataArrayInt([0,5])))
pass
+ def testIntersect2DMeshWith1DLine1(self):
+ """ Basic test for Intersect2DMeshWith1DLine: a vertical line intersecting a square. """
+ m1c = MEDCouplingCMesh()
+ coordX = DataArrayDouble([-1., 1., 2])
+ m1c.setCoordsAt(0,coordX)
+ coordY = DataArrayDouble([0., 2.])
+ m1c.setCoordsAt(1,coordY);
+ m1 = m1c.buildUnstructured()
+
+ # A simple line:
+ m2 = MEDCouplingUMesh("bla", 1)
+ coord2 = DataArrayDouble([0.,-1.0, 0.,1., 0.,3., 0.5,2.2], 4, 2)
+ conn2 = DataArrayInt([NORM_SEG2,0,1,NORM_SEG3,1,2,3])
+ connI2 = DataArrayInt([0,3,7])
+ m2.setCoords(coord2)
+ m2.setConnectivity(conn2, connI2)
+
+ # End of construction of input meshes m1bis and m2 -> start of specific part of the test
+ m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10)
+ m3.mergeNodes(1.0e-8)
+
+ self.assertEqual(3,m3.getNumberOfCells())
+ self.assertEqual(20,m3.getNumberOfNodes())
+ self.assertEqual(2,m3.getSpaceDimension())
+ # Mapping
+ exp1, exp2, expI2 = [0,0,1], [0,1,0,1], [0,2,4,4]
+ self.assertEqual(3, map1.getNumberOfTuples())
+ self.assertEqual(4, map2.getNumberOfTuples())
+ self.assertEqual(4, mapI2.getNumberOfTuples())
+ self.assertEqual(exp1, map1.getValues())
+ self.assertEqual(exp2, map2.getValues())
+ self.assertEqual(expI2, mapI2.getValues())
+
+ expConn = [32, 1, 10, 7, 11, 4, 12, 13, 14, 15, 16, 32, 10, 0, 3, 11, 7, 17, 18, 19, 14, 13, 5, 2, 1, 4, 5] # 27tpl
+ expConnI = [0, 11, 22, 27]
+ self.assertEqual(27, m3.getNodalConnectivity().getNumberOfTuples())
+ self.assertEqual(4, m3.getNodalConnectivityIndex().getNumberOfTuples())
+ self.assertEqual(expConn, m3.getNodalConnectivity().getValues())
+ self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues())
+
+ def testSwig2Intersect2DMeshWith1DLine2(self):
+ """ Star pattern (a triangle intersecting another one upside down) """
+ coords1 = DataArrayDouble([-2.,1., 2.,1., 0.,-2.], 3,2)
+ coords2 = DataArrayDouble([0.,2., 2.,-1., -2.,-1., 0.,3.], 4,2)
+ m1 = MEDCouplingUMesh("triangle", 2)
+ m2 = MEDCouplingUMesh("tri_line", 1)
+ m1.setCoords(coords1)
+ m2.setCoords(coords2)
+ m1.setConnectivity(DataArrayInt([NORM_TRI3, 0,1,2]), DataArrayInt([0,4]))
+ m2.setConnectivity(DataArrayInt([NORM_SEG2,0,1,NORM_SEG2,1,2,NORM_SEG2,2,3]), DataArrayInt([0,3,6,9]))
+
+ # End of construction of input meshes m1bis and m2 -> start of specific part of the test
+ m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10)
+ # Mapping
+ exp1, exp2, expI2 = [0,0,0,0], [2, 0,1,2, 0, 1], [0,1,4,5,6]
+ self.assertEqual(exp1, map1.getValues())
+ self.assertEqual(exp2, map2.getValues())
+ self.assertEqual(expI2, mapI2.getValues())
+ expConn = [5, 0, 8, 12, 5, 8, 7, 9, 10, 11, 12, 5, 7, 1, 9, 5, 10, 2, 11]
+ expConnI = [0, 4, 11, 15, 19]
+ self.assertEqual(expConn, m3.getNodalConnectivity().getValues())
+ self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues())
+
+ def testSwig2Intersect2DMeshWith1DLine3(self):
+ """ Line pieces ending (or fully located) in the middle of a cell """
+ m1c = MEDCouplingCMesh()
+ m1c.setCoordsAt(0,DataArrayDouble([-1., 1.]))
+ m1c.setCoordsAt(1,DataArrayDouble([-1., 1.]));
+ m1 = m1c.buildUnstructured()
+ coords2 = DataArrayDouble([0.,0., 0.,1.5, -1.5,0., 0.6,1.1, 1.1,0.6, 0.5,0.0, 0.0,-0.5], 7,2)
+ m2 = MEDCouplingUMesh("piecewise_line", 1)
+ m2.setCoords(coords2)
+ c = DataArrayInt([NORM_SEG2,0,1, NORM_SEG2,1,2, NORM_SEG2,3,4, NORM_SEG2,5,6])
+ cI = DataArrayInt([0,3,6,9,12])
+ m2.setConnectivity(c, cI)
+
+ # End of construction of input meshes m1bis and m2 -> start of specific part of the test
+ m3, map1, map2, mapI2 = MEDCouplingUMesh.Intersect2DMeshWith1DLine(m1, m2, 1e-10)
+ self.assertEqual(3,m3.getNumberOfCells())
+ # Mapping
+ exp1, exp2, expI2 = [0,0,0], [1,2, 1, 2], [0,2,3,4]
+ self.assertEqual(exp1, map1.getValues())
+ self.assertEqual(exp2, map2.getValues())
+ self.assertEqual(expI2, mapI2.getValues())
+ expConn = [5, 1, 0, 11, 13, 12, 14, 15, 5, 11, 2, 13, 5, 14, 3, 15]
+ expConnI = [0, 8, 12, 16]
+ self.assertEqual(expConn, m3.getNodalConnectivity().getValues())
+ self.assertEqual(expConnI, m3.getNodalConnectivityIndex().getValues())
+
def setUp(self):
pass
pass
if __name__ == '__main__':
unittest.main()
+
return ret;
}
+ static PyObject *Intersect2DMeshWith1DLine(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps) throw(INTERP_KERNEL::Exception)
+ {
+ DataArrayInt *cellNb1=0,*cellNb2=0, *cellNbI2=0;
+ MEDCouplingUMesh *mret=MEDCouplingUMesh::Intersect2DMeshWith1DLine(m1,m2,eps,cellNb1,cellNb2, cellNbI2);
+ PyObject *ret=PyTuple_New(4);
+ 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 ));
+ PyTuple_SetItem(ret,3,SWIG_NewPointerObj(SWIG_as_voidptr(cellNbI2),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+ return ret;
+ }
+
PyObject *buildSlice3D(PyObject *origin, PyObject *vec, double eps) const throw(INTERP_KERNEL::Exception)
{
int spaceDim=self->getSpaceDimension();
