Salome HOME
Merge from V6_main (dev of Anthony GEAY) 11/06/2013
[tools/medcoupling.git] / src / INTERP_KERNEL / Geometric2D / InterpKernelGeo2DComposedEdge.cxx
1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20
21 #include "InterpKernelGeo2DComposedEdge.hxx"
22 #include "InterpKernelGeo2DElementaryEdge.hxx"
23 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
24 #include "InterpKernelGeo2DEdgeInfLin.hxx"
25 #include "InterpKernelException.hxx"
26
27 #include <algorithm>
28 #include <memory>
29 #include <iterator>
30 #include <set>
31
32 using namespace INTERP_KERNEL;
33
34 ComposedEdge::ComposedEdge(const ComposedEdge& other)
35 {
36   for(std::list<ElementaryEdge *>::const_iterator iter=other._sub_edges.begin();iter!=other._sub_edges.end();iter++)
37     _sub_edges.push_back((*iter)->clone());
38 }
39
40 ComposedEdge::~ComposedEdge()
41 {
42   clearAll(_sub_edges.begin());
43 }
44
45 void ComposedEdge::setValueAt(int i, Edge *e, bool direction)
46 {
47   std::list<ElementaryEdge*>::iterator it=_sub_edges.begin();
48   for(int j=0;j<i;j++)
49     it++;
50   delete *it;
51   *it=new ElementaryEdge(e,direction);
52 }
53
54 struct AbsEdgeCmp
55 {
56   AbsEdgeCmp(ElementaryEdge *b):_b1(b) { }
57   bool operator()(ElementaryEdge *a) { return a->getPtr()==_b1->getPtr();}
58
59   ElementaryEdge *_b1;
60 };
61
62 double ComposedEdge::getCommonLengthWith(const ComposedEdge& other) const
63 {
64   double ret=0.;
65   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
66     {
67       if(find_if(other._sub_edges.begin(),other._sub_edges.end(),AbsEdgeCmp(*iter))!=other._sub_edges.end())
68         {
69           const ElementaryEdge *tmp=static_cast<const ElementaryEdge *>(*iter);
70           ret+=tmp->getCurveLength();
71         }
72     }
73   return ret;
74 }
75
76 void ComposedEdge::clear()
77 {
78   clearAll(_sub_edges.begin());
79   _sub_edges.clear();
80 }
81
82 void ComposedEdge::pushBack(Edge *edge, bool direction)
83 {
84   _sub_edges.push_back(new ElementaryEdge(edge,direction));
85 }
86
87 void ComposedEdge::pushBack(ElementaryEdge *elem)
88 {
89   _sub_edges.push_back(elem);
90 }
91
92 void ComposedEdge::pushBack(ComposedEdge *elem)
93 {
94   std::list<ElementaryEdge *> *elemsOfElem=elem->getListBehind();
95   _sub_edges.insert(_sub_edges.end(),elemsOfElem->begin(),elemsOfElem->end());
96 }
97
98 ElementaryEdge *ComposedEdge::operator[](int i) const
99 {
100   std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();
101   for(int ii=0;ii<i;ii++)
102     iter++;
103   return *iter;
104 }
105
106 void ComposedEdge::reverse()
107 {
108   _sub_edges.reverse();
109   for(std::list<ElementaryEdge *>::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
110     (*iter)->reverse();
111 }
112
113 bool ComposedEdge::presenceOfOn() const
114 {
115   bool ret=false;
116   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end() && !ret;iter++)
117     ret=((*iter)->getLoc()==FULL_ON_1);
118   return ret;
119 }
120
121 bool ComposedEdge::presenceOfQuadraticEdge() const
122 {
123   bool ret=false;
124   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end() && !ret;iter++)
125     {
126       Edge *e=(*iter)->getPtr();
127       if(e)
128         ret=dynamic_cast<EdgeArcCircle*>(e)!=0;
129     }
130   return ret;
131 }
132
133 void ComposedEdge::initLocations() const
134 {
135   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
136     (*iter)->initLocations();
137 }
138
139 void ComposedEdge::initLocationsWithOther(const ComposedEdge& other) const
140 {
141   std::set<Edge *> s1,s2;
142   for(std::list<ElementaryEdge *>::const_iterator it1=_sub_edges.begin();it1!=_sub_edges.end();it1++)
143     s1.insert((*it1)->getPtr());
144   for(std::list<ElementaryEdge *>::const_iterator it2=other._sub_edges.begin();it2!=other._sub_edges.end();it2++)
145     s2.insert((*it2)->getPtr());
146   initLocations();
147   other.initLocations();
148   std::vector<Edge *> s3;
149   std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::back_insert_iterator< std::vector<Edge *> >(s3));
150   for(std::vector<Edge *>::const_iterator it3=s3.begin();it3!=s3.end();it3++)
151     (*it3)->declareOn();
152 }
153
154 ComposedEdge *ComposedEdge::clone() const
155 {
156   return new ComposedEdge(*this);
157 }
158
159 bool ComposedEdge::isNodeIn(Node *n) const
160 {
161   bool ret=false;
162   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end() && !ret;iter++)
163     ret=(*iter)->isNodeIn(n);
164   return ret;
165 }
166
167 /*!
168  * This method computes the area of 'this'.
169  * By definition :
170  * \f[
171  * Area=\int_{Polygon} dS
172  * \f]
173  * Thanks to Green's theorem we have.
174  * \f[
175  * \int_{Polygon} x \cdot dS=\sum_{0 \leq i < nb of edges} -\int_{Edge_{i}}ydx=\sum_{0 \leq i < nb of edges} AreaOfZone_{Edge_{i}}
176  * \f]
177  * Where \f$ AreaOfZone_{i} \f$ is computed virtually by INTERP_KERNEL::Edge::getAreaOfZone with following formula :
178  * \f[
179  * AreaOfZone_{i}=\int_{Edge_{i}} -ydx
180  * \f]
181  */
182 double ComposedEdge::getArea() const
183 {
184   double ret=0.;
185   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
186     ret+=(*iter)->getAreaOfZone();
187   return ret;
188 }
189
190 double ComposedEdge::getPerimeter() const
191 {
192   double ret=0.;
193   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
194     ret+=(*iter)->getCurveLength();
195   return ret;
196 }
197
198 double ComposedEdge::getHydraulicDiameter() const
199 {
200   return 4*fabs(getArea())/getPerimeter();
201 }
202
203 /*!
204  * This method computes barycenter of 'this' by returning xG in bary[0] and yG in bary[1].
205  * By definition :
206  * \f[
207  * Area \cdot x_{G}=\int_{Polygon} x \cdot dS
208  * \f]
209  * \f[
210  * Area \cdot y_{G}=\int_{Polygon} y \cdot dS
211  * \f]
212  * Thanks to Green's theorem we have.
213  * \f[
214  * \int_{Polygon} x \cdot dS=\sum_{0 \leq i < nb of edges} -\int_{Edge_{i}}yxdx
215  * \f]
216  * \f[
217  * \int_{Polygon} y \cdot dS=\sum_{0 \leq i < nb of edges} -\int_{Edge_{i}}\frac{y^{2}}{2}dx
218  * \f]
219  * Area is computed using the same principle than described in INTERP_KERNEL::ComposedEdge::getArea method.
220  * \f$ -\int_{Edge_{i}}yxdx \f$ and \f$ -\int_{Edge_{i}}\frac{y^{2}}{2}dx \f$ are computed virtually with INTERP_KERNEL::Edge::getBarycenterOfZone.
221  */
222 void ComposedEdge::getBarycenter(double *bary) const
223 {
224   bary[0]=0.;
225   bary[1]=0.;
226   double area=0.;
227   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
228     {
229       (*iter)->getBarycenterOfZone(bary);
230       area+=(*iter)->getAreaOfZone();
231     }
232   bary[0]/=area;
233   bary[1]/=area;
234 }
235
236 /*!
237  * Idem ComposedEdge::getBarycenter except that the special case where _sub_edges==1 is dealt here.
238  */
239 void ComposedEdge::getBarycenterGeneral(double *bary) const
240 {
241   if(_sub_edges.empty())
242     throw INTERP_KERNEL::Exception("ComposedEdge::getBarycenterGeneral called on an empty polygon !");
243   if(_sub_edges.size()>2)
244     return getBarycenter(bary);
245   double w;
246   _sub_edges.back()->getBarycenter(bary,w);
247 }
248
249 double ComposedEdge::normalize(ComposedEdge *other, double& xBary, double& yBary)
250 {
251   Bounds b;
252   b.prepareForAggregation();
253   fillBounds(b); 
254   other->fillBounds(b);
255   double dimChar=b.getCaracteristicDim();
256   b.getBarycenter(xBary,yBary);
257   applyGlobalSimilarity(xBary,yBary,dimChar);
258   other->applyGlobalSimilarity(xBary,yBary,dimChar);
259   return dimChar;
260 }
261
262 /*!
263  * This method operates the opposite operation than ComposedEdge::applyGlobalSimilarity.
264  */
265 void ComposedEdge::unApplyGlobalSimilarityExt(ComposedEdge& other, double xBary, double yBary, double fact)
266 {
267   std::set<Node *> allNodes;
268   getAllNodes(allNodes);
269   other.getAllNodes(allNodes);
270   for(std::set<Node *>::iterator iter=allNodes.begin();iter!=allNodes.end();iter++)
271     (*iter)->unApplySimilarity(xBary,yBary,fact);
272   for(std::list<ElementaryEdge *>::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
273     (*iter)->unApplySimilarity(xBary,yBary,fact);
274   for(std::list<ElementaryEdge *>::iterator iter=other._sub_edges.begin();iter!=other._sub_edges.end();iter++)
275     (*iter)->unApplySimilarity(xBary,yBary,fact);
276 }
277
278 double ComposedEdge::normalizeExt(ComposedEdge *other, double& xBary, double& yBary)
279 {
280   Bounds b;
281   b.prepareForAggregation();
282   fillBounds(b); 
283   other->fillBounds(b);
284   double dimChar=b.getCaracteristicDim();
285   b.getBarycenter(xBary,yBary);
286   applyGlobalSimilarity2(other,xBary,yBary,dimChar);
287   return dimChar;
288 }
289
290 void ComposedEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const
291 {
292   stream.precision(10);
293   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
294     (*iter)->dumpInXfigFile(stream,resolution,box);
295 }
296
297 Node *ComposedEdge::getEndNode() const
298 {
299   return _sub_edges.back()->getEndNode();
300 }
301
302 Node *ComposedEdge::getStartNode() const
303 {
304   return _sub_edges.front()->getStartNode();
305 }
306
307 bool ComposedEdge::changeEndNodeWith(Node *node) const
308 {
309   return _sub_edges.back()->changeEndNodeWith(node);
310 }
311
312 bool ComposedEdge::changeStartNodeWith(Node *node) const
313 {
314   return _sub_edges.front()->changeStartNodeWith(node);
315 }
316
317 void ComposedEdge::fillBounds(Bounds& output) const
318 {
319   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
320     (*iter)->fillBounds(output);
321 }
322
323 /*!
324  * \b WARNING : applies similarity \b ONLY on edges without any change on Nodes. To perform a global similarity call applyGlobalSimilarity.
325  */
326 void ComposedEdge::applySimilarity(double xBary, double yBary, double dimChar)
327 {
328   for(std::list<ElementaryEdge *>::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
329     (*iter)->applySimilarity(xBary,yBary,dimChar);
330 }
331
332 /*!
333  * Perform Similarity transformation on all elements of this Nodes and Edges.
334  */
335 void ComposedEdge::applyGlobalSimilarity(double xBary, double yBary, double dimChar)
336 {
337   std::set<Node *> allNodes;
338   getAllNodes(allNodes);
339   for(std::set<Node *>::iterator iter=allNodes.begin();iter!=allNodes.end();iter++)
340     (*iter)->applySimilarity(xBary,yBary,dimChar);
341   for(std::list<ElementaryEdge *>::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
342     (*iter)->applySimilarity(xBary,yBary,dimChar);
343 }
344
345 /*!
346  * Perform Similarity transformation on all elements of this Nodes and Edges on 'this' and 'other'.
347  * Nodes can be shared between 'this' and 'other'.
348  */
349 void ComposedEdge::applyGlobalSimilarity2(ComposedEdge *other, double xBary, double yBary, double dimChar)
350 {
351   std::set<Node *> allNodes;
352   getAllNodes(allNodes);
353   std::set<Node *> allNodes2;
354   other->getAllNodes(allNodes2);
355   for(std::set<Node *>::const_iterator it=allNodes2.begin();it!=allNodes2.end();it++)
356     if(allNodes.find(*it)!=allNodes.end())
357       (*it)->declareOn();
358   allNodes.insert(allNodes2.begin(),allNodes2.end());
359   for(std::set<Node *>::iterator iter=allNodes.begin();iter!=allNodes.end();iter++)
360     (*iter)->applySimilarity(xBary,yBary,dimChar);
361   for(std::list<ElementaryEdge *>::iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
362     (*iter)->applySimilarity(xBary,yBary,dimChar);
363   for(std::list<ElementaryEdge *>::iterator iter=other->_sub_edges.begin();iter!=other->_sub_edges.end();iter++)
364     (*iter)->applySimilarity(xBary,yBary,dimChar);
365 }
366
367 /*!
368  * This method append to param 'partConsidered' the part of length of subedges IN or ON.
369  * @param partConsidered INOUT param.
370  */
371 void ComposedEdge::dispatchPerimeter(double& partConsidered) const
372 {
373   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
374     {
375       TypeOfEdgeLocInPolygon loc=(*iter)->getLoc();
376       if(loc==FULL_IN_1 || loc==FULL_ON_1)
377         partConsidered+=(*iter)->getCurveLength();
378     }
379 }
380
381 /*!
382  * Idem dispatchPerimeterExcl except that when a subedge is declared as ON this subedge is counted in commonPart.
383  */
384 void ComposedEdge::dispatchPerimeterExcl(double& partConsidered, double& commonPart) const
385 {
386   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
387     {
388       TypeOfEdgeLocInPolygon loc=(*iter)->getLoc();
389       if(loc==FULL_IN_1)
390         partConsidered+=(*iter)->getCurveLength();
391       if(loc==FULL_ON_1)
392         commonPart+=(*iter)->getCurveLength();
393     }
394 }
395
396 void ComposedEdge::getAllNodes(std::set<Node *>& output) const
397 {
398   std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();
399   for(;iter!=_sub_edges.end();iter++)
400     (*iter)->getAllNodes(output);
401 }
402
403 void ComposedEdge::getBarycenter(double *bary, double& weigh) const
404 {
405   weigh=0.; bary[0]=0.; bary[1]=0.;
406   double tmp1,tmp2[2];
407   for(std::list<ElementaryEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
408     {
409       (*iter)->getBarycenter(tmp2,tmp1);
410       weigh+=tmp1;
411       bary[0]+=tmp1*tmp2[0];
412       bary[1]+=tmp1*tmp2[1];
413     }
414   bary[0]/=weigh;
415   bary[1]/=weigh;
416 }
417
418 bool ComposedEdge::isInOrOut(Node *nodeToTest) const
419 {
420   Bounds b; b.prepareForAggregation();
421   fillBounds(b);
422   if(b.nearlyWhere((*nodeToTest)[0],(*nodeToTest)[1])==OUT)
423     return false;
424   // searching for e1
425   std::set<Node *> nodes;
426   getAllNodes(nodes);
427   std::set<double> radialDistributionOfNodes;
428   std::set<Node *>::const_iterator iter;
429   for(iter=nodes.begin();iter!=nodes.end();iter++)
430     radialDistributionOfNodes.insert(nodeToTest->getSlope(*(*iter)));
431   std::vector<double> radialDistrib(radialDistributionOfNodes.begin(),radialDistributionOfNodes.end());
432   radialDistributionOfNodes.clear();
433   std::vector<double> radialDistrib2(radialDistrib.size());
434   copy(radialDistrib.begin()+1,radialDistrib.end(),radialDistrib2.begin());
435   radialDistrib2.back()=M_PI+radialDistrib.front();
436   std::vector<double> radialDistrib3(radialDistrib.size());
437   std::transform(radialDistrib2.begin(),radialDistrib2.end(),radialDistrib.begin(),radialDistrib3.begin(),std::minus<double>());
438   std::vector<double>::iterator iter3=max_element(radialDistrib3.begin(),radialDistrib3.end());
439   int i=iter3-radialDistrib3.begin();
440   // ok for e1 - Let's go.
441   EdgeInfLin *e1=new EdgeInfLin(nodeToTest,radialDistrib[i]+radialDistrib3[i]/2.);
442   double ref=e1->getCharactValue(*nodeToTest);
443   std::set< IntersectElement > inOutSwitch;
444   for(std::list<ElementaryEdge *>::const_iterator iter4=_sub_edges.begin();iter4!=_sub_edges.end();iter4++)
445     {
446       ElementaryEdge *val=(*iter4);
447       if(val)
448         {
449           Edge *e=val->getPtr();
450           std::auto_ptr<EdgeIntersector> intersc(Edge::BuildIntersectorWith(e1,e));
451           bool obviousNoIntersection,areOverlapped;
452           intersc->areOverlappedOrOnlyColinears(0,obviousNoIntersection,areOverlapped);
453           if(obviousNoIntersection)
454             {
455               continue;
456             }
457           if(!areOverlapped)
458             {
459               std::list< IntersectElement > listOfIntesc=intersc->getIntersectionsCharacteristicVal();
460               for(std::list< IntersectElement >::iterator iter2=listOfIntesc.begin();iter2!=listOfIntesc.end();iter2++)
461                 if((*iter2).isIncludedByBoth())
462                   inOutSwitch.insert(*iter2);
463               }
464           //if overlapped we can forget
465         }
466       else
467         throw Exception("Invalid use of ComposedEdge::isInOrOut : only one level supported !");
468     }
469   e1->decrRef();
470   bool ret=false;
471   for(std::set< IntersectElement >::iterator iter4=inOutSwitch.begin();iter4!=inOutSwitch.end();iter4++)
472     {
473       if((*iter4).getVal1()<ref)
474         {
475           if((*iter4).getNodeOnly()->getLoc()==ON_1)
476             ret=!ret;
477         }
478       else
479         break;
480     }
481   return ret;
482 }
483
484 /*bool ComposedEdge::isInOrOut(Node *aNodeOn, Node *nodeToTest) const
485 {
486   
487   EdgeInfLin *e1=new EdgeInfLin(aNodeOn,nodeToTest);
488   double ref=e1->getCharactValue(*nodeToTest);
489   set< IntersectElement > inOutSwitch;
490   for(vector<AbstractEdge *>::const_iterator iter=_sub_edges.begin();iter!=_sub_edges.end();iter++)
491     {
492       ElementaryEdge *val=dynamic_cast<ElementaryEdge *>(*iter);
493       if(val)
494         {
495           Edge *e=val->getPtr();
496           auto_ptr<Intersector> intersc(Edge::buildIntersectorWith(e1,e));
497           bool obviousNoIntersection,areOverlapped;
498           intersc->areOverlappedOrOnlyColinears(0,obviousNoIntersection,areOverlapped);
499           if(obviousNoIntersection)
500             {
501               continue;
502             }
503           if(!areOverlapped)
504             {
505               list< IntersectElement > listOfIntesc=intersc->getIntersectionsCharacteristicVal();
506               for(list< IntersectElement >::iterator iter2=listOfIntesc.begin();iter2!=listOfIntesc.end();iter2++)
507                 if((*iter2).isIncludedByBoth())
508                   inOutSwitch.insert(*iter2);
509               }
510           //if overlapped we can forget
511         }
512       else
513         throw Exception("Invalid use of ComposedEdge::isInOrOut : only one level supported !");
514     }
515   e1->decrRef();
516   bool ret=false;
517   for(set< IntersectElement >::iterator iter=inOutSwitch.begin();iter!=inOutSwitch.end();iter++)
518     {
519       if((*iter).getVal1()<ref)
520         {
521           if((*iter).getNodeOnly()->getLoc()==ON_1)
522             ret=!ret;
523         }
524       else
525         break;
526     }
527   return ret;
528 }*/
529
530 bool ComposedEdge::getDirection() const
531 {
532   throw Exception("ComposedEdge::getDirection : no sense");
533 }
534
535 bool ComposedEdge::intresincEqCoarse(const Edge *other) const
536 {
537   if(_sub_edges.size()!=1)
538     return false;
539   return _sub_edges.front()->intresincEqCoarse(other);
540 }
541
542 void ComposedEdge::clearAll(std::list<ElementaryEdge *>::iterator startToDel)
543 {
544   for(std::list<ElementaryEdge *>::iterator iter=startToDel;iter!=_sub_edges.end();iter++)
545     delete (*iter);
546 }