src/VTKViewer/VTKViewer_ConvexTool.cxx

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