-// Copyright (C) 2007-2012 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2013 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
+// Author : Anthony Geay (CEA/DEN)
#include "InterpKernelGeo2DQuadraticPolygon.hxx"
#include "InterpKernelGeo2DElementaryEdge.hxx"
#include "NormalizedUnstructuredMesh.hxx"
#include <fstream>
+#include <sstream>
#include <iomanip>
#include <cstring>
#include <limits>
*/
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,
- const std::vector< std::vector<int> >& colinear1)
+ const std::vector< std::vector<int> >& colinear1,
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
{
std::size_t nbOfSeg=std::distance(descBg,descEnd);
for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
{
bool direct=descBg[i]>0;
int edgeId=abs(descBg[i])-1;//current edge id of pol2
+ std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
+ if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end())
+ {
+ bool sameDir=(it1!=alreadyExistingIn2.end());
+ const std::vector<INTERP_KERNEL::ElementaryEdge *>& edgesAlreadyBuilt=sameDir?(*it1).second:(*it2).second;
+ if(sameDir)
+ {
+ 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()));
+ }
+ }
+ else
+ {
+ for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_reverse_iterator it4=edgesAlreadyBuilt.rbegin();it4!=edgesAlreadyBuilt.rend();it4++)
+ {
+ Edge *ee=(*it4)->getPtr(); ee->incrRef();
+ pushBack(new ElementaryEdge(ee,!(*it4)->getDirection()));
+ }
+ }
+ continue;
+ }
bool directos=colinear1[edgeId].empty();
std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
int offset1=0;
}
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();
appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
+ std::size_t newSz=_sub_edges.size();
+ std::size_t zeSz=newSz-oldSz;
+ alreadyExistingIn2[descBg[i]].resize(zeSz);
+ std::list<ElementaryEdge *>::const_reverse_iterator it5=_sub_edges.rbegin();
+ for(std::size_t p=0;p<zeSz;p++,it5++)
+ alreadyExistingIn2[descBg[i]][zeSz-p-1]=*it5;
}
else
{//there is subpart of edge 'edgeId' of pol2 inside pol1
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();it++)
+ 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;
Node *end=(*mapp.find(idEnd)).second;
ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
pushBack(e);
+ alreadyExistingIn2[descBg[i]].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(); ee->declareOn();
- pushBack(new ElementaryEdge(ee,!(direct1^direction11)));
+ ee->incrRef();
+ ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
+ pushBack(e2);
+ alreadyExistingIn2[descBg[i]].push_back(e2);
}
}
}
}
}
+/*!
+ * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
+ */
+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 nbOfSeg=std::distance(descBg,descEnd);
+ for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
+ {
+ bool direct=descBg[i]>0;
+ int edgeId=abs(descBg[i])-1;//current edge id of pol2
+ const std::vector<int>& c=colinear1[edgeId];
+ if(c.empty())
+ continue;
+ 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++)
+ {
+ 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=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
+ int idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
+ 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
+ }
+ }
+}
+
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
{
int nbOfNodesInPg=0;
/*!
* 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'.
+ * @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
*/
-void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, 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'
- other.performLocatingOperation(*this);
+ other.performLocatingOperationSlow(*this);
std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
{
(*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,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);
+ }
+ }
nbThis.push_back(idThis);
nbOther.push_back(idOther);
delete *it;
}
}
+void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const
+{
+ IteratorOnComposedEdge it(&pol2);
+ for(it.first();!it.finished();it.next())
+ {
+ ElementaryEdge *cur=it.current();
+ cur->locateFullyMySelfAbsolute(*this);
+ }
+}
+
/*!
* Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned.
*
std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
{
std::list<QuadraticPolygon *> ret;
- IteratorOnComposedEdge it((ComposedEdge *)this);
+ IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(this));
int nbOfTurns=recursiveSize();
int i=0;
if(!it.goToNextInOn(false,i,nbOfTurns))
return iter;
return iEnd;
}
+
+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)
+{
+ pol1.initLocations();
+ for(std::set<Edge *>::const_iterator it=notUsedInPol1.begin();it!=notUsedInPol1.end();it++)
+ { (*it)->initLocs(); (*it)->declareOn(); }
+ for(std::set<Edge *>::const_iterator it=edgesInPol2OnBoundary.begin();it!=edgesInPol2OnBoundary.end();it++)
+ { (*it)->initLocs(); (*it)->declareIn(); }
+ ////
+ std::set<Edge *> notUsedInPol1L(notUsedInPol1);
+ IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1));
+ int sz=pol1.size();
+ 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);
+ return ;
+ }
+ while(!notUsedInPol1L.empty())
+ {
+ for(int i=0;i<sz && (it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1);i++)
+ it.nextLoop();
+ if(it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1)
+ throw INTERP_KERNEL::Exception("Presence of a target polygon fully included in source polygon ! The partition of this leads to a non simply connex cell (with hole) ! Impossible ! Such resulting cell cannot be stored in MED cell format !");
+ QuadraticPolygon *tmp1=new QuadraticPolygon;
+ do
+ {
+ Edge *ee=it.current()->getPtr();
+ if(ee->getLoc()==FULL_ON_1)
+ {
+ ee->incrRef(); notUsedInPol1L.erase(ee);
+ tmp1->pushBack(new ElementaryEdge(ee,it.current()->getDirection()));
+ }
+ it.nextLoop();
+ }
+ while(it.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
+ pol1Zip.push_back(tmp1);
+ }
+ ////
+ std::list<QuadraticPolygon *> retPolsUnderContruction;
+ std::list<Edge *> edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end());
+ std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed;
+ std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0;
+ for(std::list<QuadraticPolygon *>::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++)
+ nbOfTurn+=(*itMNT)->size();
+ maxNbOfTurn=maxNbOfTurn*nbOfTurn; maxNbOfTurn*=maxNbOfTurn;
+ nbOfTurn=0;
+ while(nbOfTurn<maxNbOfTurn && ((!pol1Zip.empty() || !edgesInPol2OnBoundaryL.empty())))
+ {
+ for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();)
+ {
+ if((*it1)->getStartNode()==(*it1)->getEndNode())
+ {
+ it1++;
+ continue;
+ }
+ Node *curN=(*it1)->getEndNode();
+ bool smthHappened=false;
+ for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
+ {
+ if(curN==(*it2)->getStartNode())
+ { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
+ else if(curN==(*it2)->getEndNode())
+ { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
+ else
+ it2++;
+ }
+ if(smthHappened)
+ {
+ for(std::list<QuadraticPolygon *>::iterator it3=pol1Zip.begin();it3!=pol1Zip.end();)
+ {
+ if(curN==(*it3)->getStartNode())
+ {
+ for(std::list<ElementaryEdge *>::const_iterator it4=(*it3)->_sub_edges.begin();it4!=(*it3)->_sub_edges.end();it4++)
+ { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*it1)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
+ smthHappened=true;
+ pol1ZipConsumed[*it1].push_back(*it3);
+ curN=(*it3)->getEndNode();
+ it3=pol1Zip.erase(it3);
+ }
+ else
+ it3++;
+ }
+ }
+ if(!smthHappened)
+ {
+ for(std::list<ElementaryEdge *>::const_iterator it5=(*it1)->_sub_edges.begin();it5!=(*it1)->_sub_edges.end();it5++)
+ {
+ Edge *ee=(*it5)->getPtr();
+ if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end())
+ edgesInPol2OnBoundaryL.push_back(ee);
+ }
+ for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
+ pol1Zip.push_front(*it6);
+ pol1ZipConsumed.erase(*it1);
+ delete *it1;
+ it1=retPolsUnderContruction.erase(it1);
+ }
+ }
+ if(!pol1Zip.empty())
+ {
+ QuadraticPolygon *tmp=new QuadraticPolygon;
+ QuadraticPolygon *first=*(pol1Zip.begin());
+ for(std::list<ElementaryEdge *>::const_iterator it4=first->_sub_edges.begin();it4!=first->_sub_edges.end();it4++)
+ { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); tmp->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
+ pol1ZipConsumed[tmp].push_back(first);
+ retPolsUnderContruction.push_back(tmp);
+ pol1Zip.erase(pol1Zip.begin());
+ }
+ nbOfTurn++;
+ }
+ if(nbOfTurn==maxNbOfTurn)
+ {
+ std::ostringstream oss; oss << "Error during reconstruction of residual of cell ! It appears that either source or/and target mesh is/are not conform !";
+ oss << " Number of turns is = " << nbOfTurn << " !";
+ throw INTERP_KERNEL::Exception(oss.str().c_str());
+ }
+ for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();it1++)
+ {
+ if((*it1)->getStartNode()==(*it1)->getEndNode())
+ {
+ (*it1)->appendCrudeData(mapp,0.,0.,1.,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;
+ }
+ }
+}
