src/VTKViewer/VTKViewer_ConvexTool.cxx

   1 // Copyright (C) 2005  OPEN CASCADE, CEA/DEN, EDF R&D, PRINCIPIA 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/
  18 //
  19 //  This library is distributed in the hope that it will be useful,
  20 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  21 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  22 //  Lesser General Public License for more details.
  23 //
  24 //  You should have received a copy of the GNU Lesser General Public
  25 //  License along with this library; if not, write to the Free Software
  26 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  27 //
  28 //  See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
  29
  30 #include "VTKViewer_ConvexTool.h"
  31
  32 #include <vtkUnstructuredGrid.h>
  33 #include <vtkTriangle.h>
  34 #include <vtkConvexPointSet.h>
  35 #include <vtkMath.h>
  36
  37 #include <set>
  38 #include <iterator>
  39 #include <algorithm>
  40 #include <math.h>
  41
  42 typedef std::set<vtkIdType> TUIDS; // unique ids
  43 typedef std::map<vtkIdType,TUIDS> TPTOIDS; // id points -> unique ids
  44
  45 namespace CONVEX_TOOL
  46 {
  47
  48 static float FACE_ANGLE_TOLERANCE=1.5;
  49
  50 #ifdef _DEBUG_
  51   static int MYDEBUG = 0;
  52 #else
  53   static int MYDEBUG = 0;
  54 #endif
  55
  56 /*! \fn static void GetCenter(vtkUnstructuredGrid* theGrid,TCell theptIds,float *center)
  57  * \brief Calculation of geometry center.
  58  * \param theGrid - vtkUnstructuredGrid cell.
  59  * \param theptIds - point ids.
  60  * \retval center - output array[3] with coordinates of geometry center.
  61  */
  62 static void GetCenter(vtkUnstructuredGrid* theGrid,TCell theptIds,float center[3])
  63 {
  64   float p[3];
  65   center[0] = center[1] = center[2] = 0.0;
  66
  67   int numIds = theptIds.size();
  68   if (numIds == 0) return;
  69
  70   // get the center of the cell
  71   for (int i = 0; i < numIds; i++)
  72   {
  73     theGrid->GetPoint(theptIds[i], p);
  74     for (int j = 0; j < 3; j++)
  75     {
  76       center[j] += p[j];
  77     }
  78   }
  79   for (int j = 0; j < 3; j++)
  80   {
  81     center[j] /= numIds;
  82   }
  83 }
  84
  85 /*! \fn static void ReverseIds(TCell &theIds)
  86  * \brief Reverse ids.
  87  * \param theIds - points ids.
  88  * \retval theIds - example input:(1,2,3,4) -> output:(4,3,2,1)
  89  */
  90 static void ReverseIds(TCell &theIds)
  91 {
  92   int i;
  93   vtkIdType tmp;
  94   vtkIdType npts=theIds.size();
  95
  96   for(i=0;i<(npts/2);i++){
  97     tmp = theIds[i];
  98     theIds[i] = theIds[npts-i-1];
  99     theIds[npts-i-1] = tmp;
 100   }
 101 }
 102
 103 /*! \fn void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output)
 104  * \brief Caclulation of connected faces (faceId -> (faceId1,faceId2, ...))
 105  * \param p2faces - point to faces ids map.
 106  * \param f2points - faces to points ids map.
 107  * \retval face2face_output - faces to faces ids map.
 108  */
 109 void GetFriends(const TPTOIDS p2faces,const TCellArray f2points,TPTOIDS& face2face_output)
 110 {
 111   TCellArray::const_iterator f2pIter = f2points.begin();
 112
 113   for( ; f2pIter!=f2points.end() ; f2pIter++ ){
 114     vtkIdType faceId = f2pIter->first;
 115     TCell face_points = f2pIter->second;
 116     int nb_face_points = face_points.size();
 117
 118     vtkIdType id1;
 119     vtkIdType id2;
 120     TPTOIDS::const_iterator faces1;
 121     TPTOIDS::const_iterator faces2;
 122
 123     id1 = face_points[0];
 124     faces1 = p2faces.find(id1);
 125
 126     TUIDS output_faces;
 127
 128     for(int i=1 ; i<nb_face_points ; i++ ){
 129
 130       id2 = face_points[i];
 131
 132       faces2 = p2faces.find(id2);
 133
 134       std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
 135         std::inserter(output_faces,output_faces.begin()));
 136
 137       id1 = id2;
 138       faces1 = faces2;
 139     }
 140     id1 = face_points[0];
 141     faces1 = p2faces.find(id1);
 142     std::set_intersection(faces1->second.begin(), faces1->second.end(), faces2->second.begin(), faces2->second.end(),
 143       std::inserter(output_faces,output_faces.begin()));
 144
 145     output_faces.erase(faceId); // erase the face id for which we found friends
 146
 147     if(MYDEBUG){
 148       cout << "fId[" << faceId <<"]: ";
 149       std::copy(output_faces.begin(), output_faces.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 150       cout << endl;
 151     }
 152
 153     face2face_output[faceId] = output_faces;
 154   }
 155 }
 156
 157 /*! \fn bool IsConnectedFacesOnOnePlane( vtkUnstructuredGrid* theGrid,vtkIdType theFId1, vtkIdType theFId2,TUIDS FpIds1, TUIDS FpIds2 )
 158  * \brief Check is connected faces on one plane.
 159  * \param theGrid - vtkUnstructuredGrid
 160  * \param theFId1 - id of first face
 161  * \param theFId2 - id of second face
 162  * \param FpIds1  - first face points ids.
 163  * \param FpIds2  - second face points ids.
 164  * \return TRUE if two faces on one plane, else FALSE.
 165  */
 166 bool IsConnectedFacesOnOnePlane( vtkUnstructuredGrid* theGrid,
 167                                  vtkIdType theFId1, vtkIdType theFId2,
 168                                  TUIDS FpIds1, TUIDS FpIds2 )
 169 {
 170   bool status = false;
 171   TUIDS common_ids;
 172   std::set_intersection(FpIds1.begin(), FpIds1.end(), FpIds2.begin(), FpIds2.end(),
 173                         std::inserter(common_ids,common_ids.begin()));
 174
 175   /*           Number of common ids = 2 (A1,A2)
 176
 177
 178                 _ _ _ _ _      _ _ _ _    vectors:
 179                |         \   /         |   v1 {A2,A1}
 180                           \ /              v2 {A1,B1}
 181                |           | A2        |   v3 {A1,C1}
 182                            |
 183                |           |           |
 184                            |
 185                |           | A1        |
 186                           / \
 187                |_ _ _ _ _/   \_ _ _ _ _|
 188                B2        B1   C1        C2
 189
 190    */
 191   TUIDS::iterator common_iter = common_ids.begin();
 192   if(common_ids.size() == 2){
 193     TUIDS::iterator loc_id1_0 = FpIds1.find(*(common_iter));
 194     common_iter++;
 195     TUIDS::iterator loc_id1_1 = FpIds1.find(*(common_iter));
 196
 197     TUIDS::iterator loc_id2_0 = FpIds1.begin();
 198     TUIDS::iterator loc_id2_1 = FpIds2.begin();
 199
 200     vtkIdType A1 = *loc_id1_0;
 201     vtkIdType A2 = *loc_id1_1;
 202     vtkIdType B1;
 203     vtkIdType C1;
 204
 205     for(;loc_id2_0!=FpIds1.end();loc_id2_0++)
 206       if(*loc_id2_0 != A1 && *loc_id2_0!= A2){
 207         B1 = *loc_id2_0;
 208         break;
 209       }
 210     for(;loc_id2_1!=FpIds2.end();loc_id2_1++)
 211       if(*loc_id2_1 != A1 && *loc_id2_1!= A2){
 212         C1 = *loc_id2_1;
 213         break;
 214       }
 215     if(MYDEBUG) cout <<endl;
 216     if(MYDEBUG) cout << "FId_1="<<theFId1<<" FId_2="<<theFId2<<endl;
 217     if(MYDEBUG) cout << "   A1="<<A1<<" A2="<<A2<<" B1="<<B1<<" C1="<<C1<<" ->";
 218     float *p[4];
 219     float v1[3],v2[3],v3[3];
 220     p[0] = theGrid->GetPoint(A1);
 221     p[1] = theGrid->GetPoint(A2);
 222     p[2] = theGrid->GetPoint(B1);
 223     p[3] = theGrid->GetPoint(C1);
 224
 225     for(int i=0;i<3;i++){
 226       v1[i] = p[1][i] - p[0][i];
 227       v2[i] = p[2][i] - p[0][i];
 228       v3[i] = p[3][i] - p[0][i];
 229     }
 230
 231
 232     float vec_b1[3];
 233     vtkMath::Cross(v2,v1,vec_b1);
 234     float vec_b2[3];
 235     vtkMath::Cross(v1,v3,vec_b2);
 236
 237     float b1 = vtkMath::Norm(vec_b1);
 238
 239     float b2 = vtkMath::Norm(vec_b2);
 240     float aCos = vtkMath::Dot(vec_b1,vec_b2)/(b1*b2);
 241
 242     float angle=90.0;
 243     angle = aCos>=1.0 ? 0.0 : 180*acosf(aCos)/vtkMath::Pi();
 244
 245     if( angle <= FACE_ANGLE_TOLERANCE)
 246       status = true;
 247     if (MYDEBUG){
 248       for(int k=0;k<4;k++){
 249         cout << " (";
 250         for(int j=0;j<2;j++){
 251           cout << p[k][j] << ",";
 252         }
 253         cout << p[k][2] << ")   ";
 254       }
 255       cout << "angle="<<angle<<" status="<<status<<endl;
 256     }
 257
 258   } else if (common_ids.size() >2){
 259     // not implemented yet
 260     if(MYDEBUG) cout << "Warning! VTKViewer_ConvexTool::IsConnectedFacesOnOnePlane()";
 261   } else {
 262     // one or no connection ... continue
 263   }
 264
 265   return status;
 266 }
 267
 268 /*! \fn void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,TUIDS &new_faces, TCell &new_faces_v2 )
 269  * \brief Calculate faces which on one plane.
 270  * \param theFaces -
 271  * \param faceId -
 272  * \param new_faces -
 273  * \param new_faces_v2 -
 274  */
 275 void GetAllFacesOnOnePlane( TPTOIDS theFaces, vtkIdType faceId,
 276                             TUIDS &new_faces, TCell &new_faces_v2 )
 277 {
 278   if (new_faces.find(faceId) != new_faces.end()) return;
 279
 280   new_faces.insert(new_faces.begin(),faceId);
 281   new_faces_v2.push_back(faceId);
 282
 283   TPTOIDS::const_iterator aIter1 = theFaces.find(faceId);
 284   if(aIter1!=theFaces.end()){
 285     TUIDS::const_iterator aIter2 = (aIter1->second).begin();
 286     for(;aIter2!=(aIter1->second).end();aIter2++){
 287       if (new_faces.find(*aIter2) != new_faces.end()) continue;
 288       GetAllFacesOnOnePlane(theFaces,*aIter2,
 289                             new_faces,new_faces_v2); // recurvise
 290     }
 291   }
 292   return;
 293 }
 294
 295 /*! \fn void GetSumm(TCell v1,TCell v2,TCell &output)
 296  * \brief Gluing two faces (gluing points ids)
 297  * \param v1 - first face
 298  * \param v2 - second face
 299  * \param output - output face.
 300  */
 301 void GetSumm(TCell v1,TCell v2,TCell &output)
 302 {
 303   output.clear();
 304
 305   if(MYDEBUG) cout << "========================================="<<endl;
 306   if(MYDEBUG) cout << "v1:";
 307   if(MYDEBUG) std::copy(v1.begin(), v1.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 308   if(MYDEBUG) cout << "\tv2:";
 309   if(MYDEBUG) std::copy(v2.begin(), v2.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 310   if(MYDEBUG) cout << endl;
 311
 312   TUIDS v1_set;
 313   std::copy(v1.begin(), v1.end(), std::inserter(v1_set,v1_set.begin()));
 314   TUIDS v2_set;
 315   std::copy(v2.begin(), v2.end(), std::inserter(v2_set,v2_set.begin()));
 316   TUIDS tmpIntersection;
 317   std::set_intersection(v1_set.begin(),v1_set.end(),v2_set.begin(),v2_set.end(), std::inserter(tmpIntersection,tmpIntersection.begin()));
 318   if(MYDEBUG) std::copy(tmpIntersection.begin(),tmpIntersection.end(), std::ostream_iterator<vtkIdType>(cout, " "));
 319   if(MYDEBUG) cout << endl;
 320
 321   if(tmpIntersection.size() < 2)
 322     if(MYDEBUG) cout << __FILE__ << "[" << __LINE__ << "]: Warning ! Wrong ids" << endl;
 323
 324   TCell::iterator v1_iter = v1.begin();
 325
 326   for(;v1_iter!=v1.end();v1_iter++){
 327
 328     vtkIdType curr_id = *v1_iter;
 329
 330     output.push_back(curr_id);
 331
 332     if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 333       TCell::iterator v1_iter_tmp;
 334       v1_iter_tmp = v1_iter;
 335       v1_iter++;
 336
 337       if(v1_iter==v1.end()) v1_iter=v1.begin();
 338
 339       curr_id = *v1_iter;
 340
 341       if(tmpIntersection.find(curr_id) != tmpIntersection.end()){
 342         TCell::iterator v2_iter = v2.begin();
 343         for(;v2_iter!=v2.end();v2_iter++){
 344           vtkIdType v2_id = *v2_iter;
 345           if(tmpIntersection.find(v2_id) == tmpIntersection.end())
 346             output.push_back(v2_id);
 347         }
 348       }
 349
 350       v1_iter = v1_iter_tmp;
 351       curr_id = *v1_iter;
 352
 353     }
 354   }
 355
 356   if(MYDEBUG) cout << "Result: " ;
 357   if(MYDEBUG) std::copy(output.begin(),output.end(),std::ostream_iterator<vtkIdType>(cout, " "));
 358   if(MYDEBUG) cout << endl;
 359 }
 360
 361 void GetPolygonalFaces(vtkUnstructuredGrid* theGrid,int cellId,TCellArray &outputCellArray)
 362 {
 363   if (theGrid->GetCellType(cellId) == VTK_CONVEX_POINT_SET){
 364     // get vtkCell type = VTK_CONVEX_POINT_SET
 365     if(vtkConvexPointSet* convex = static_cast<vtkConvexPointSet*>(theGrid->GetCell(cellId))){
 366       TCellArray f2points;
 367       float convex_center[3]; // convex center point coorinat
 368       int aNbFaces = convex->GetNumberOfFaces();
 369       int numPts = convex->GetNumberOfPoints();
 370       TCell convex_ids;
 371       TPTOIDS p2faces; // key=pointId , value facesIds set
 372
 373       for(int i=0;i<numPts;i++)
 374           convex_ids.push_back(convex->GetPointId(i));
 375
 376       GetCenter(theGrid,convex_ids,convex_center);
 377
 378       for (vtkIdType faceId=0; faceId < aNbFaces; faceId++){
 379         vtkCell *aFace = convex->GetFace(faceId);
 380         int numFacePts = aFace->GetNumberOfPoints();
 381         TCell aIds;
 382
 383         int i = 0;
 384         for(i=0;i<numFacePts;i++)
 385           aIds.push_back(aFace->GetPointId(i));
 386
 387         float v_a[3],v_b[3],v_convex2face[3]; // vectors
 388         float *id_0,*id_1,*id_n;
 389         /*=============================================
 390                         ,+- - - -  _
 391                    _   / id_n   |  v_b {id_0,id_n}
 392                   v_b /            _
 393                      /          |  v_a {id_0,id_1}
 394                     /
 395                    /            |
 396                   + id_0
 397                    \
 398                   _ \           |
 399                  v_a \
 400                       \ id_1    |
 401                        "+- - - -
 402
 403          ============================================*/
 404         id_0 = theGrid->GetPoint(aIds[0]);
 405         id_1 = theGrid->GetPoint(aIds[1]);
 406         id_n = theGrid->GetPoint(aIds[numFacePts-1]);
 407
 408         for(i=0;i<3;i++){
 409           v_a[i] = id_1[i] - id_0[i];
 410           v_b[i] = id_n[i] - id_0[i];
 411           v_convex2face[i] = id_0[i] - convex_center[i];
 412         }
 413
 414         if (vtkMath::Determinant3x3(v_a,v_b,v_convex2face) < 0){
 415           ReverseIds(aIds);
 416         }
 417
 418         for(i=0;i<(int)aIds.size();i++){
 419           TUIDS &acell = p2faces[aIds[i]];
 420           acell.insert(faceId);
 421         }
 422
 423         f2points[faceId] = aIds;
 424
 425       }
 426
 427       TPTOIDS face2face;
 428       GetFriends(p2faces,f2points,face2face);
 429
 430       TPTOIDS face2points;
 431
 432       // copy TCellArray::f2points to TPTOIDS::face2points
 433       for(TCellArray::iterator f2points_iter=f2points.begin();
 434           f2points_iter!=f2points.end();
 435           f2points_iter++){
 436
 437         TUIDS tmp;
 438         for(TCell::iterator points_iter=(f2points_iter->second).begin();
 439             points_iter!=(f2points_iter->second).end();
 440             points_iter++)
 441           tmp.insert(*points_iter);
 442
 443         face2points[f2points_iter->first] = tmp;
 444       } // end copy
 445
 446
 447       TPTOIDS new_face2faces; // which connected and in one plane
 448
 449       TPTOIDS::const_iterator aF2FIter = face2face.begin();
 450       for(;aF2FIter!=face2face.end();aF2FIter++){
 451         vtkIdType f_key = aF2FIter->first;
 452         TUIDS &faces = new_face2faces[f_key];
 453         //faces.insert(f_key);
 454         TUIDS f_friends = aF2FIter->second;
 455         TUIDS::const_iterator a_friends_iter = f_friends.begin();
 456         for(;a_friends_iter!=f_friends.end();a_friends_iter++){
 457           vtkIdType friend_id = *a_friends_iter;
 458           if( IsConnectedFacesOnOnePlane(theGrid,f_key,friend_id,
 459                                         (face2points.find(f_key))->second,
 460                                         (face2points.find(friend_id))->second)){
 461             faces.insert(friend_id);
 462           } // end if
 463
 464         } // end a_friends_iter
 465       } // end aF2FIter
 466
 467       if(MYDEBUG)
 468       {
 469         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 470         cout << "Connected faces and on plane:" << endl;
 471         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 472           cout << "Group ["<<new_face2face_iter->first<<"] :";
 473           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 474           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 475             cout << " " << *new_faces_iter ;
 476           cout << endl;
 477         }
 478       }
 479
 480       TPTOIDS output_newid2face;
 481       TCellArray output_newid2face_v2;
 482       {
 483         TUIDS already_in;
 484         TUIDS already_in_tmp;
 485         int k=0;
 486         TPTOIDS::const_iterator new_face2face_iter = new_face2faces.begin();
 487         for(;new_face2face_iter!=new_face2faces.end();new_face2face_iter++){
 488           if(already_in.find(new_face2face_iter->first) != already_in.end())
 489             continue;
 490           if(new_face2face_iter->second.size() > 1)
 491             continue;
 492
 493           TCell &tmp_v2 = output_newid2face_v2[k];
 494           tmp_v2.push_back(new_face2face_iter->first);
 495           already_in.insert(new_face2face_iter->first);
 496
 497           TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 498           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 499             if(already_in.find(*new_faces_iter) != already_in.end()) continue;
 500             already_in.insert(*new_faces_iter);
 501
 502             already_in_tmp.clear();
 503             already_in_tmp.insert(new_face2face_iter->first);
 504
 505             TUIDS &tmp = output_newid2face[k];
 506             GetAllFacesOnOnePlane(new_face2faces,*new_faces_iter,
 507                                   already_in_tmp,tmp_v2);
 508
 509             for(TUIDS::const_iterator aIter=already_in_tmp.begin();
 510                 aIter!=already_in_tmp.end();
 511                 aIter++)
 512               {
 513                 already_in.insert(*aIter);
 514                 tmp.insert(*aIter);
 515               }
 516           }
 517           k++;
 518         }
 519       }
 520
 521       if(MYDEBUG) {
 522         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 523         {
 524           TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 525           for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 526             cout << "Group ["<<new_face2face_iter->first<<"] :";
 527             TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 528             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 529               cout << " " << *new_faces_iter ;
 530             cout << endl;
 531           }
 532         }
 533         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 534         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 535         cout << "+++++++++++++++++++++++ ++ ++ ++++++++++++++++++++++++++++"<<endl;
 536         cout << "+++++++++++++++++++++++++   ++++++++++++++++++++++++++++++"<<endl;
 537         cout << "++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++"<<endl;
 538         {
 539           TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 540           for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 541             cout << "Group ["<<new_face2face_iter->first<<"] :";
 542             TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 543             for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++)
 544               cout << " " << *new_faces_iter ;
 545             cout << endl;
 546           }
 547         }
 548       }
 549       TCellArray output_new_face2ids;
 550 //       {
 551 //      TPTOIDS::const_iterator new_face2face_iter = output_newid2face.begin();
 552 //      for(;new_face2face_iter!=output_newid2face.end();new_face2face_iter++){
 553
 554 //        vtkIdType new_faceId = new_face2face_iter->first;
 555 //        TUIDS::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 556 //        vtkIdType fId0 = *new_faces_iter;
 557 //        TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 558 //        TCell pIds0 = pIds0_iter->second;
 559 //        TCell &output = output_new_face2ids[new_faceId];
 560 //        new_faces_iter++;
 561 //        if(new_face2face_iter->second.size() > 2 ){}
 562 //        for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 563
 564 //          vtkIdType faceId = *new_faces_iter;
 565 //          // find how much nodes in face (f2points)
 566 //          TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 567 //          TCell pIds = pIds_iter->second;
 568
 569 //          GetSumm(pIds0,pIds,output);
 570 //          pIds0 = output;
 571
 572 //        } // end new_faces_iter
 573
 574 //      } // new_face2face_iter
 575 //       }
 576
 577       {
 578         TCellArray::const_iterator new_face2face_iter = output_newid2face_v2.begin();
 579         for(;new_face2face_iter!=output_newid2face_v2.end();new_face2face_iter++){
 580
 581           vtkIdType new_faceId = new_face2face_iter->first;
 582           TCell::const_iterator new_faces_iter = (new_face2face_iter->second).begin();
 583           vtkIdType fId0 = *new_faces_iter;
 584           TCellArray::const_iterator pIds0_iter = f2points.find(fId0);
 585           TCell pIds0 = pIds0_iter->second;
 586           TCell &output = output_new_face2ids[new_faceId];
 587           new_faces_iter++;
 588           if(new_face2face_iter->second.size() == 1 ){
 589             TCellArray::const_iterator pIds_iter = f2points.find(fId0);
 590             TCell pIds = pIds_iter->second;
 591             output = pIds;
 592             continue;
 593           }
 594           for(;new_faces_iter!=(new_face2face_iter->second).end();new_faces_iter++){
 595
 596             vtkIdType faceId = *new_faces_iter;
 597             // find how much nodes in face (f2points)
 598             TCellArray::const_iterator pIds_iter = f2points.find(faceId);
 599             TCell pIds = pIds_iter->second;
 600
 601             GetSumm(pIds0,pIds,output);
 602             pIds0 = output;
 603
 604           } // end new_faces_iter
 605
 606         } // new_face2face_iter
 607       }
 608
 609       if(MYDEBUG) {
 610         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 611         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 612         cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 613       }
 614       outputCellArray = output_new_face2ids;//f2points;
 615     }
 616   } else {
 617     // not implemented
 618   }
 619 }
 620 }