Save Triangle Mesh in std::vector using SMPtools

Hello,

I`m trying to save my 3 million triangle polydata with 2 scalar fields in a std::vector with SMPtools. But I get an error in the
vtkDataArrayTupleRange_AOS. I tryed to print the begin and end tuples to check what was happening.

Thank you in advance,

Rafael Scatena

        std::vector<Triangle> h_triangles(numTriangles);
          
        ///////    SPM
    
          std::cerr << "PRISMAMESH: No of Triangles: " << trianglePolyData->GetNumberOfCells()<< std::endl;
      
     
          vtkSMPTools::For(0, numTriangles, [&](vtkIdType begin, vtkIdType end)       {




      

          std::cerr << "beginTuple: " << begin << ", endTuple: " << end << std::endl;


          for (; begin < end; ++begin) {
          
              vtkCell* cell = trianglePolyData->GetCell(begin);
              vtkTriangle* triangle = vtkTriangle::SafeDownCast(cell);



              vtkIdType pointId0, pointId1, pointId2; 
          
              // Extract indices
              pointId0 = triangle->GetPointId(0);
              pointId1 = triangle->GetPointId(1);
              pointId2 = triangle->GetPointId(2);




              h_triangles[begin].triangle_id = static_cast<int64_t>(begin);
                



      		vtkIdType pointIds[3] = {pointId0, pointId1, pointId2};  double coordinates[3];
              
              points->GetPoint(pointIds[0], coordinates);
              
              h_triangles[begin].p1x = static_cast<double>(coordinates[0]);
              h_triangles[begin].p1y = static_cast<double>(coordinates[1]);
              h_triangles[begin].p1z = static_cast<double>(coordinates[2]);
              

              points->GetPoint(pointIds[1], coordinates);

              h_triangles[begin].p2x = static_cast<double>(coordinates[0]);
              h_triangles[begin].p2y = static_cast<double>(coordinates[1]);
              h_triangles[begin].p2z = static_cast<double>(coordinates[2]);

              points->GetPoint(pointIds[2], coordinates);

              h_triangles[begin].p3x = static_cast<double>(coordinates[0]);
              h_triangles[begin].p3y = static_cast<double>(coordinates[1]);
              h_triangles[begin].p3z = static_cast<double>(coordinates[2]);


              h_triangles[begin].Material = static_cast<int>(materialArray->GetTuple1(begin));
              h_triangles[begin].Body     = static_cast<int>(bodiesArray->GetTuple1(begin));            
              
                              
              }});
              
              std::cout << "Triangle List Ready ! "<< std::endl;
                              
              return h_triangles;};

The error issued is this one:

: /home/3484681/miniconda3/envs/vtkenv/vtk-master/source2/VTK/Common/Core/vtkDataArrayTupleRange_AOS.h:921: vtk::detail::TupleRange<vtkAOSDataArrayTemplate<ValueType>, TupleSize>::TupleRange(vtk::detail::TupleRange<vtkAOSDataArrayTemplate<ValueType>, TupleSize>::ArrayType*, vtk::TupleIdType, vtk::TupleIdType) [with ValueType = double; int TupleSize = 3; vtk::detail::TupleRange<vtkAOSDataArrayTemplate<ValueType>, TupleSize>::ArrayType = vtkAOSDataArrayTemplate<double>; vtk::TupleIdType = long long int]: Assertion `endTuple >= 0 && endTuple <= this->Array->GetNumberOfTuples()' failed.
, endTuple: 648349
61165, endTuple: 73398beginTuple: 770679, endTuple: 782912
758446Abortado (imagem do núcleo gravada)

This line is not thread-safe, and is probably the reason for the error:

The reason it isn’t thread-safe is that a vtkPolyData object only has one vtkTriangle object, and when you call this method, it updates that triangle’s IDs according to the requested cell Id. So the various threads end up fighting over what the triangle’s IDs should be.


For efficiency, a good way to get the triangles is to get the cell array:

vtkCellArray* triangles = trianglePolyData->GetPolys();

Then you can index into the cell array using this method:

void vtkCellArray::GetCellAtId(vtkIdType cellId, vtkIdList* pts);

This requires a thread-local vtkIdList object to hold the point IDs.

Thank you

vtkSmartPointer<vtkCellArray> polys = trianglePolyData->GetPolys();

          vtkSmartPointer<vtkIntArray> materialArray = vtkIntArray::SafeDownCast(trianglePolyData->GetCellData()->GetArray("MATERIAL"));
          vtkSmartPointer<vtkIntArray> bodiesArray = vtkIntArray::SafeDownCast(trianglePolyData->GetCellData()->GetArray("IBODY"));

          
          if (!materialArray || !bodiesArray) {std::cerr << "Error: Missing scalar data." << std::endl;return {};}

          
          int numTriangles = trianglePolyData->GetNumberOfCells();
       
       
          std::vector<Triangle> h_triangles(numTriangles);
          
          ///////    SPM
      
          std::cerr << "PRISMAMESH: No of Triangles: " << trianglePolyData->GetNumberOfCells()<< std::endl;
      
     
          vtkSMPTools::For(0, numTriangles-1, [&](vtkIdType begin, vtkIdType end){


          for (vtkIdType begin=0; begin < end ; ++begin) {
          

              vtkSmartPointer<vtkIdList> pointIdList = vtkSmartPointer<vtkIdList>::New(); // Create vtkIdList

               polys->GetCellAtId(begin, pointIdList);
              
              int i = static_cast<int>(begin);

              // Extract the points and store them in your custom structure
              vtkIdType pointId0, pointId1, pointId2; 
          
              // Extract indices
              pointId0 = pointIdList->GetId(0);
              pointId1 = pointIdList->GetId(1);
              pointId2 = pointIdList->GetId(2);

              std::cerr << "PRISMAMESH: Cell Id:  "<<i << std::endl;

              h_triangles[i].triangle_id = static_cast<int64_t>(i);


      		vtkIdType pointIds[3];  double coordinates[3];
              
              points->GetPoint(pointId0, coordinates);
              
             
              h_triangles[i].p1x = static_cast<double>(coordinates[0]);
              h_triangles[i].p1y = static_cast<double>(coordinates[1]);
              h_triangles[i].p1z = static_cast<double>(coordinates[2]);
     
              points->GetPoint(pointId1, coordinates);

              h_triangles[i].p2x = static_cast<double>(coordinates[0]);
              h_triangles[i].p2y = static_cast<double>(coordinates[1]);
              h_triangles[i].p2z = static_cast<double>(coordinates[2]);

              points->GetPoint(pointId2, coordinates);

              h_triangles[i].p3x = static_cast<double>(coordinates[0]);
              h_triangles[i].p3y = static_cast<double>(coordinates[1]);
              h_triangles[i].p3z = static_cast<double>(coordinates[2]);

              h_triangles[i].Material = static_cast<int>(materialArray->GetTuple1(i));
              h_triangles[i].Body     = static_cast<int>(bodiesArray->GetTuple1(i));            
              
                              
              }
              //}std::cerr << "Invalid loop bounds: begin = " << begin << ", end = " << end << std::endl;
              return;});
              
              std::cout << "Triangle List Ready ! "<< std::endl;