<commit_before># standard python modules
import re

# extensions
import vtk
import numpy

from icqBaseShape import BaseShape
from icqUniformGrid import uniformGrid
from icqTriangulation import Triangulation

class CompositeShape(BaseShape):

  def __init__(self,):
    """
    Constructor
    """
    BaseShape.__init__(self,)

    self.expression = ''
    self.argShapes = []

  def assemble(self, expression, shapes):
    """
    Assemble shapes into a composite object
    @param expression expression of the shape objects, $0 for shape 0, $1 for shape 1, etc. 
           operations include + (union), * (intersection), and - (removal)
    """
    self.expression = expression
    self.argShapes = shapes
    # replace $1 by self.argShape[0]
    for i in range(len(self.argShapes)):
      self.expression = re.sub('\${0}'.format(i), 
        'self.argShapes[{0}].evaluate(pts)'.format(i), self.expression)

  def evaluate(self, pts):
    """
    Evaluate the inside/outside function 
    @param pts array of points
    """
    return eval(self.expression)

  def computeSurfaceMeshes(self, maxTriArea):
    """
    Compute the surface meshes
    @param maxTriArea maximum triangle area
    """
  
    # compute cell normals
    self.computeSurfaceNormals()

    #
    # find all the points inside the volume
    #

    # move point out by a small amount to avoid floating point comparision
    # issues
    eps = 1.e-6
    ONE_THIRD = 1./3.

    print '*** self.surfaceMeshes = ', self.surfaceMeshes
    numFaces = len(self.surfaceMeshes)
    insidePointList = []
    totalNumPoints = 0
    loBound = numpy.array([float('inf')] * 3)
    hiBound = numpy.array([-float('inf')] * 3)
    for shp in self.argShapes:
      numFaces = len(shp.surfaceMeshes)
      for iFace in range(numFaces):
        print '*** iFace = ', iFace
        face = shp.surfaceMeshes[iFace]
        shp.computeSurfaceNormals()
        normals = shp.surfaceNormals[iFace]
        print '*** normals = ', normals
        print '*** face = ', face
        print '*** points = ', shp.points
        verts = shp.points[face] 
        displVerts = verts.copy()
        displVerts[:, 0, :] += eps*normals
        displVerts[:, 1, :] += eps*normals
        displVerts[:, 2, :] += eps*normals
        inPts = verts[self.evaluate(displVerts) > 0]
        totalNumPoints += inPts.shape[0]
        insidePointList.append(inPts)
        xs, ys, zs = inPts[:, 0], inPts[:, 1], inPts[:, 2]
        print '*** xs, ys, zs = ', xs, ys, zs
        loBound[0] = min(loBound[0], xs.min())
        loBound[1] = min(loBound[1], ys.min())
        loBound[2] = min(loBound[2], zs.min())
        hiBound[0] = max(hiBound[0], xs.max())
        hiBound[1] = max(hiBound[1], ys.max())
        hiBound[2] = max(hiBound[2], zs.max())

    # add more points by sampling the volume containing the boundary surfaces
    cellLength = numpy.sqrt(2*maxTriArea)
    ns = numpy.array([ int( (hiBound[i] - loBound[i])/cellLength + 0.5) for i in range(3) ])
    xx, yy, zz = uniformGrid(loBound, hiBound, ns)
    npts = len(xx.flat)
    gridPts = numpy.zeros( (npts, 3), numpy.float64 )
    gridPts[:, 0] = xx.flat
    gridPts[:, 1] = yy.flat
    gridPts[:, 2] = zz.flat

    # only add the points inside
    gridPts = gridPts[ self.evaluate(gridPts) > 0 ]
    totalNumPoints += gridPts.shape[0]
    insidePointList.append(gridPts)

    # create a single list of points
    insidePoints = numpy.zeros( (totalNumPoints, 3), numpy.float64 )
    iBeg = 0
    for pts in insidePointList:
      iEnd = iBeg + pts.shape[0]
      insidePoints[iBeg:iEnd, :] = pts
      iBeg = iEnd 

    # tessellate
    tri = Triangulation()
    tri.setInputPoints(insidePoints)
    tri.triangulate()

    ugrid = tri.delny.GetOutput()

    # remove the cells that are outside the closed object
    pts = ugrid.GetPoints()
    cells = ugrid.GetCells()
    ptIds = vtk.vtkIdList()
    print '*** cells = ', cells
    numCells = cells.GetNumberOfCells()
    volumeMesh = numpy.zeros( (numCells, 4), numpy.int )
    centroidPoints = numpy.zeros( (numCells, 3,), numpy.float64 )
    for i in range(numCells):
      cell = cells.GetCell(i, ptIds)
      i0, i1, i2, i3 = ptIds.GetId(0), ptIds.GetId(1), ptIds.GetId(2), ptIds.GetId(3)
      p0, p1, p2, p3 = pts.GetPoint(i0), pts.GetPoint(i1), pts.GetPoint(i2), pts.GetPoint(i3)
      centroidPoints[i, :] += numpy.array(p0)
      centroidPoints[i, :] += numpy.array(p1)
      centroidPoints[i, :] += numpy.array(p2)
      centroidPoints[i, :] += numpy.array(p3)
      centroidPoints[i, :] *= 0.25
      volumeMesh[i, :] = i0, i1, i2, i3

    validCells = self.evaluate(centroidPoints)
    print '*** validCells = ', validCells

    ugrid2 = vtk.vtkUnstructuredGrid()
    ugrid2.SetPoints(pts)
    numCells2 = int(validCells.sum())
    ugrid2.Allocate(numCells2, 1)
    for i in range(numCells):
      if validCells[i]:
        cell = cells.GetCell(i, ptIds)
        ugrid2.InsertNextCell(vtk.VTK_TETRA, 4, volumeMesh[i, :])

    # apply filter to extract boundary cell faces

    print '*** ugrid2 = ', ugrid2

    surf = vtk.vtkGeometryFilter()
    if vtk.VTK_MAJOR_VERSION >= 6:
      surf.SetInputData(ugrid2)
    else:
      surf.SetInput(ugrid2)
    surf.Update()

    print '*** surf = ', surf

    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(surf.GetOutputPort())

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)

    ren = vtk.vtkRenderer()
    win = vtk.vtkRenderWindow()
    iren = vtk.vtkRenderWindowInteractor()

    ren.AddActor(actor)
    win.AddRenderer(ren)
    iren.SetRenderWindow(win)
    win.SetSize(500, 500)<commit_msg> GetNextCell works better<commit_after># standard python modules
import re

# extensions
import vtk
import numpy

from icqBaseShape import BaseShape
from icqUniformGrid import uniformGrid
from icqTriangulation import Triangulation

class CompositeShape(BaseShape):

  def __init__(self,):
    """
    Constructor
    """
    BaseShape.__init__(self,)

    self.expression = ''
    self.argShapes = []

  def assemble(self, expression, shapes):
    """
    Assemble shapes into a composite object
    @param expression expression of the shape objects, $0 for shape 0, $1 for shape 1, etc. 
           operations include + (union), * (intersection), and - (removal)
    """
    self.expression = expression
    self.argShapes = shapes
    # replace $1 by self.argShape[0]
    for i in range(len(self.argShapes)):
      self.expression = re.sub('\${0}'.format(i), 
        'self.argShapes[{0}].evaluate(pts)'.format(i), self.expression)

  def evaluate(self, pts):
    """
    Evaluate the inside/outside function 
    @param pts array of points
    """
    return eval(self.expression)

  def computeSurfaceMeshes(self, maxTriArea):
    """
    Compute the surface meshes
    @param maxTriArea maximum triangle area
    """
  
    # compute cell normals
    self.computeSurfaceNormals()

    #
    # find all the points inside the volume
    #

    # move point out by a small amount to avoid floating point comparision
    # issues
    eps = 1.e-6
    ONE_THIRD = 1./3.

    print '*** self.surfaceMeshes = ', self.surfaceMeshes
    numFaces = len(self.surfaceMeshes)
    insidePointList = []
    totalNumPoints = 0
    loBound = numpy.array([float('inf')] * 3)
    hiBound = numpy.array([-float('inf')] * 3)
    for shp in self.argShapes:
      numFaces = len(shp.surfaceMeshes)
      for iFace in range(numFaces):
        print '*** iFace = ', iFace
        face = shp.surfaceMeshes[iFace]
        shp.computeSurfaceNormals()
        normals = shp.surfaceNormals[iFace]
        print '*** normals = ', normals
        print '*** face = ', face
        print '*** points = ', shp.points
        verts = shp.points[face] 
        displVerts = verts.copy()
        displVerts[:, 0, :] += eps*normals
        displVerts[:, 1, :] += eps*normals
        displVerts[:, 2, :] += eps*normals
        inPts = verts[self.evaluate(displVerts) > 0]
        totalNumPoints += inPts.shape[0]
        insidePointList.append(inPts)
        xs, ys, zs = inPts[:, 0], inPts[:, 1], inPts[:, 2]
        print '*** xs, ys, zs = ', xs, ys, zs
        loBound[0] = min(loBound[0], xs.min())
        loBound[1] = min(loBound[1], ys.min())
        loBound[2] = min(loBound[2], zs.min())
        hiBound[0] = max(hiBound[0], xs.max())
        hiBound[1] = max(hiBound[1], ys.max())
        hiBound[2] = max(hiBound[2], zs.max())

    # add more points by sampling the volume containing the boundary surfaces
    cellLength = numpy.sqrt(2*maxTriArea)
    ns = numpy.array([ int( (hiBound[i] - loBound[i])/cellLength + 0.5) for i in range(3) ])
    xx, yy, zz = uniformGrid(loBound, hiBound, ns)
    npts = len(xx.flat)
    gridPts = numpy.zeros( (npts, 3), numpy.float64 )
    gridPts[:, 0] = xx.flat
    gridPts[:, 1] = yy.flat
    gridPts[:, 2] = zz.flat

    # only add the points inside
    gridPts = gridPts[ self.evaluate(gridPts) > 0 ]
    totalNumPoints += gridPts.shape[0]
    insidePointList.append(gridPts)

    # create a single list of points
    insidePoints = numpy.zeros( (totalNumPoints, 3), numpy.float64 )
    iBeg = 0
    for pts in insidePointList:
      iEnd = iBeg + pts.shape[0]
      insidePoints[iBeg:iEnd, :] = pts
      iBeg = iEnd 

    # tessellate
    tri = Triangulation()
    tri.setInputPoints(insidePoints)
    tri.triangulate()

    ugrid = tri.delny.GetOutput()

    # remove the cells that are outside the closed object
    pts = ugrid.GetPoints()
    cells = ugrid.GetCells()
    ptIds = vtk.vtkIdList()
    print '*** cells = ', cells
    numCells = cells.GetNumberOfCells()
    volumeMesh = numpy.zeros( (numCells, 4), numpy.int )
    centroidPoints = numpy.zeros( (numCells, 3,), numpy.float64 )
    for i in range(numCells):
      cell = cells.GetCell(i, ptIds)
      i0, i1, i2, i3 = ptIds.GetId(0), ptIds.GetId(1), ptIds.GetId(2), ptIds.GetId(3)
      p0, p1, p2, p3 = pts.GetPoint(i0), pts.GetPoint(i1), pts.GetPoint(i2), pts.GetPoint(i3)
      centroidPoints[i, :] += numpy.array(p0)
      centroidPoints[i, :] += numpy.array(p1)
      centroidPoints[i, :] += numpy.array(p2)
      centroidPoints[i, :] += numpy.array(p3)
      centroidPoints[i, :] *= 0.25
      volumeMesh[i, :] = i0, i1, i2, i3

    validCells = self.evaluate(centroidPoints)
    print '*** validCells = ', validCells

    ugrid2 = vtk.vtkUnstructuredGrid()
    ugrid2.SetPoints(pts)
    numCells2 = int(validCells.sum())
    ugrid2.Allocate(numCells2, 1)
    cells.InitTraversal()
    for i in range(numCells):
      if validCells[i]:
        cell = cells.GetNextCell(ptIds)
        print '*** number of Ids = ', ptIds.GetNumberOfIds()
        ugrid2.InsertNextCell(vtk.VTK_TETRA, ptIds)

    # apply filter to extract boundary cell faces

    print '*** ugrid = ', ugrid
    print '*** ugrid2 = ', ugrid2
    
    surf = vtk.vtkGeometryFilter()
    if vtk.VTK_MAJOR_VERSION >= 6:
      surf.SetInputData(ugrid2)
    else:
      surf.SetInput(ugrid2)
    surf.Update()

    print '*** surf = ', surf

    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(surf.GetOutputPort())

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)

    ren = vtk.vtkRenderer()
    win = vtk.vtkRenderWindow()
    iren = vtk.vtkRenderWindowInteractor()

    ren.AddActor(actor)
    win.AddRenderer(ren)
    iren.SetRenderWindow(win)
    win.SetSize(500, 500)
{'fraction_non_alphanumeric': 0.1352920478536242, 'fraction_numerical': 0.017857142857142856, 'mean_word_length': 2.8305256064690028, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '912325', 'n_tokens_mistral': 4159, 'n_tokens_neox': 3992, 'n_words': 901}