Free Form Deformation with PyGem

Marco Tezzele (https://github.com/mtezzele) and Filippo Salmoiraghi (https://github.com/fsalmoir) develop the PyGem project. They contributed the following blog post to explain their project. Many thanks !

logo_PyGeM_small

PyGeM (Python Geometrical Morphing) is a python package using Free Form Deformation and Radial Basis Functions to parametrize and morph complex geometries. It is ideally suited for actual industrial problems, since it allows to handle:

  • Computer Aided Design files (in .iges and .stl formats)
  • Mesh files (in .unv and OpenFOAM formats)
  • Output files (in .vtk format)

Since it has to integrate itself in the industrial workflow we have chosen Python. You can find the project on GitHub (https://github.com/mathLab/PyGeM).

Just to have an overview of the possibilities offered by PyGeM we show some examples:

ffd_MCY_hull
MCY hull: morphing of the exhaust gasses devices starting from an industrial .stl file.

compare_geometries_with_lattice

ffdDrivAer

DrivAer model: morphing of the bumper starting from an OpenFOAM mesh file.

FFD introduction

Here we explain in details the Free Form Deformation technique in order to better understand the industrial problem. In the documentation of the package you can also find details about the Radial Basis Functions interpolation technique.

Free Form Deformation is a technique for the efficient, smooth and accurate geometrical parametrization. It has been proposed the first time in Sederberg, Thomas W., and Scott R. Parry. “Free-form deformation of solid geometric models.” ACM SIGGRAPH computer graphics 20.4 (1986): 151-160. It has been developed for computer graphics, but in the last decades it has proven as a very useful tool for fluid dynamics shape optimization.

It consists in three different step:

ffd_scheme
Then we can compose the three maps in the following way:
composition

Then every point inside the FFD box changes it position according to
points

where b denotes the Bernstein polynomials. We improve the traditional version by allowing a rotation of the FFD lattice in order to give more flexibility to the tool.
You can try to add more shapes to the lattice to allow more and more involved transformations.

FFD on IGES: the role of PythonOCC

In order to perform FFD on iges files we exploit the PythonOCC package. We mainly use it as a tool to parse and write iges files. Every CAD file is made of several surfaces, described by NURBS functions. We perform the FFD directly on the poles (or control points). In particular we extract the poles coordinates, we move them according to the user preferences (set through a parameters file) and then we write back the new poles coordinates in order to visualize the deformed geometry.

PythonOCC is heavily used in the igeshandler module. In the following we show only the parser and the writer.

Code for the iges parsing:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
# read in the IGES file
reader = IGESControl_Reader()
reader.ReadFile(self.infile)
reader.TransferRoots()
shape = reader.Shape()
 
# cycle on the faces to get the control points
# init some quantities
n_faces = 0
control_point_position = [0]
faces_explorer = TopExp_Explorer(shape, TopAbs_FACE)
mesh_points = np.zeros(shape=(0, 3))
 
while faces_explorer.More():
    # performing some conversions to get the right format (BSplineSurface)
    iges_face = OCC.TopoDS.topods_Face(faces_explorer.Current())
    iges_nurbs_converter = BRepBuilderAPI_NurbsConvert(iges_face)
    iges_nurbs_converter.Perform(iges_face)
    nurbs_face = iges_nurbs_converter.Shape()
    brep_face = BRep_Tool.Surface(OCC.TopoDS.topods_Face(nurbs_face))
    bspline_face = geomconvert_SurfaceToBSplineSurface(brep_face)
 
    # openCascade object
    occ_face = bspline_face.GetObject()
 
    # extract the Control Points of each face
    n_poles_u = occ_face.NbUPoles()
    n_poles_v = occ_face.NbVPoles()
    control_polygon_coordinates = np.zeros(shape=(n_poles_u * n_poles_v, 3))
    
    # cycle over the poles to get their coordinates
    i = 0
    for pole_u_direction in xrange(n_poles_u):
        for pole_v_direction in xrange(n_poles_v):
            control_point_coordinates = occ_face.Pole(pole_u_direction+1, pole_v_direction+1)
            control_polygon_coordinates[i, :] = [control_point_coordinates.X(), \
            control_point_coordinates.Y(), control_point_coordinates.Z()]
            i += 1
 
# pushing the control points coordinates to the mesh_points array (used for FFD)
mesh_points = np.append(mesh_points, control_polygon_coordinates, axis=0)
control_point_position.append(control_point_position[-1] + n_poles_u * n_poles_v)
 
n_faces += 1
faces_explorer.Next()
# read in the IGES file
reader = IGESControl_Reader()
reader.ReadFile(self.infile)
reader.TransferRoots()
shape = reader.Shape()

# cycle on the faces to get the control points
# init some quantities
n_faces = 0
control_point_position = [0]
faces_explorer = TopExp_Explorer(shape, TopAbs_FACE)
mesh_points = np.zeros(shape=(0, 3))

while faces_explorer.More():
    # performing some conversions to get the right format (BSplineSurface)
    iges_face = OCC.TopoDS.topods_Face(faces_explorer.Current())
    iges_nurbs_converter = BRepBuilderAPI_NurbsConvert(iges_face)
    iges_nurbs_converter.Perform(iges_face)
    nurbs_face = iges_nurbs_converter.Shape()
    brep_face = BRep_Tool.Surface(OCC.TopoDS.topods_Face(nurbs_face))
    bspline_face = geomconvert_SurfaceToBSplineSurface(brep_face)

    # openCascade object
    occ_face = bspline_face.GetObject()

    # extract the Control Points of each face
    n_poles_u = occ_face.NbUPoles()
    n_poles_v = occ_face.NbVPoles()
    control_polygon_coordinates = np.zeros(shape=(n_poles_u * n_poles_v, 3))
    
    # cycle over the poles to get their coordinates
    i = 0
    for pole_u_direction in xrange(n_poles_u):
        for pole_v_direction in xrange(n_poles_v):
            control_point_coordinates = occ_face.Pole(pole_u_direction+1, pole_v_direction+1)
            control_polygon_coordinates[i, :] = [control_point_coordinates.X(), \
            control_point_coordinates.Y(), control_point_coordinates.Z()]
            i += 1

# pushing the control points coordinates to the mesh_points array (used for FFD)
mesh_points = np.append(mesh_points, control_polygon_coordinates, axis=0)
control_point_position.append(control_point_position[-1] + n_poles_u * n_poles_v)

n_faces += 1
faces_explorer.Next()

And here we have the writer:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# init the ouput file writer
writer = IGESControl_Writer()
 
# read in the IGES file
reader = IGESControl_Reader()
reader.ReadFile(self.infile)
reader.TransferRoots()
shape_read = reader.Shape()
 
# cycle on the faces to update the control points position
# init some quantities
faces_explorer = TopExp_Explorer(shape_read, TopAbs_FACE)
n_faces = 0
control_point_position = self._control_point_position
 
compound_builder = BRep_Builder()
compound = OCC.TopoDS.TopoDS_Compound()
compound_builder.MakeCompound(compound)
 
while faces_explorer.More():
    # similar to the parser method
    iges_face = OCC.TopoDS.topods_Face(faces_explorer.Current())
    iges_nurbs_converter = BRepBuilderAPI_NurbsConvert(iges_face)
    iges_nurbs_converter.Perform(iges_face)
    nurbs_face = iges_nurbs_converter.Shape()
    face_aux = OCC.TopoDS.topods_Face(nurbs_face)
    brep_face = BRep_Tool.Surface(OCC.TopoDS.topods_Face(nurbs_face))
    bspline_face = geomconvert_SurfaceToBSplineSurface(brep_face)
    occ_face = bspline_face.GetObject()
    
    n_poles_u = occ_face.NbUPoles()
    n_poles_v = occ_face.NbVPoles()
 
    i = 0
    for pole_u_direction in xrange(n_poles_u):
        for pole_v_direction in xrange(n_poles_v):
            control_point_coordinates = mesh_points[i+control_point_position[n_faces], :]
            point_xyz = gp_XYZ(control_point_coordinates[0], control_point_coordinates[1], \
            control_point_coordinates[2])
            gp_point = gp_Pnt(point_xyz)
            occ_face.SetPole(pole_u_direction+1, pole_v_direction+1, gp_point)
            i += 1
 
    # construct the deformed wire for the trimmed surfaces
    wire_maker = BRepBuilderAPI_MakeWire()
    tol = ShapeFix_ShapeTolerance()
    brep = BRepBuilderAPI_MakeFace(occ_face.GetHandle(), self.tolerance).Face()
    brep_face = BRep_Tool.Surface(brep)
 
    # cycle on the edges
    edge_explorer = TopExp_Explorer(nurbs_face, TopAbs_EDGE)
    while edge_explorer.More():
        edge = OCC.TopoDS.topods_Edge(edge_explorer.Current())
        # edge in the (u,v) coordinates
        edge_uv_coordinates = OCC.BRep.BRep_Tool.CurveOnSurface(edge, face_aux)
        # evaluating the new edge: same (u,v) coordinates, but different (x,y,x) ones
        edge_phis_coordinates_aux = BRepBuilderAPI_MakeEdge(edge_uv_coordinates[0], brep_face)
        edge_phis_coordinates = edge_phis_coordinates_aux.Edge()
        tol.SetTolerance(edge_phis_coordinates, self.tolerance)
        wire_maker.Add(edge_phis_coordinates)
        edge_explorer.Next()
 
    # grouping the edges in a wire
    wire = wire_maker.Wire()
 
    # trimming the surfaces
    brep_surf = BRepBuilderAPI_MakeFace(occ_face.GetHandle(), wire).Shape()
    compound_builder.Add(compound, brep_surf)
    n_faces += 1
    faces_explorer.Next()
 
writer.AddShape(compound)
 
writer.Write(self.outfile)
# init the ouput file writer
writer = IGESControl_Writer()

# read in the IGES file
reader = IGESControl_Reader()
reader.ReadFile(self.infile)
reader.TransferRoots()
shape_read = reader.Shape()

# cycle on the faces to update the control points position
# init some quantities
faces_explorer = TopExp_Explorer(shape_read, TopAbs_FACE)
n_faces = 0
control_point_position = self._control_point_position

compound_builder = BRep_Builder()
compound = OCC.TopoDS.TopoDS_Compound()
compound_builder.MakeCompound(compound)

while faces_explorer.More():
    # similar to the parser method
    iges_face = OCC.TopoDS.topods_Face(faces_explorer.Current())
    iges_nurbs_converter = BRepBuilderAPI_NurbsConvert(iges_face)
    iges_nurbs_converter.Perform(iges_face)
    nurbs_face = iges_nurbs_converter.Shape()
    face_aux = OCC.TopoDS.topods_Face(nurbs_face)
    brep_face = BRep_Tool.Surface(OCC.TopoDS.topods_Face(nurbs_face))
    bspline_face = geomconvert_SurfaceToBSplineSurface(brep_face)
    occ_face = bspline_face.GetObject()
    
    n_poles_u = occ_face.NbUPoles()
    n_poles_v = occ_face.NbVPoles()

    i = 0
    for pole_u_direction in xrange(n_poles_u):
        for pole_v_direction in xrange(n_poles_v):
            control_point_coordinates = mesh_points[i+control_point_position[n_faces], :]
            point_xyz = gp_XYZ(control_point_coordinates[0], control_point_coordinates[1], \
            control_point_coordinates[2])
            gp_point = gp_Pnt(point_xyz)
            occ_face.SetPole(pole_u_direction+1, pole_v_direction+1, gp_point)
            i += 1

    # construct the deformed wire for the trimmed surfaces
    wire_maker = BRepBuilderAPI_MakeWire()
    tol = ShapeFix_ShapeTolerance()
    brep = BRepBuilderAPI_MakeFace(occ_face.GetHandle(), self.tolerance).Face()
    brep_face = BRep_Tool.Surface(brep)

    # cycle on the edges
    edge_explorer = TopExp_Explorer(nurbs_face, TopAbs_EDGE)
    while edge_explorer.More():
        edge = OCC.TopoDS.topods_Edge(edge_explorer.Current())
        # edge in the (u,v) coordinates
        edge_uv_coordinates = OCC.BRep.BRep_Tool.CurveOnSurface(edge, face_aux)
        # evaluating the new edge: same (u,v) coordinates, but different (x,y,x) ones
        edge_phis_coordinates_aux = BRepBuilderAPI_MakeEdge(edge_uv_coordinates[0], brep_face)
        edge_phis_coordinates = edge_phis_coordinates_aux.Edge()
        tol.SetTolerance(edge_phis_coordinates, self.tolerance)
        wire_maker.Add(edge_phis_coordinates)
        edge_explorer.Next()

    # grouping the edges in a wire
    wire = wire_maker.Wire()

    # trimming the surfaces
    brep_surf = BRepBuilderAPI_MakeFace(occ_face.GetHandle(), wire).Shape()
    compound_builder.Add(compound, brep_surf)
    n_faces += 1
    faces_explorer.Next()

writer.AddShape(compound)

writer.Write(self.outfile)

Here we have a simple example of the input/output procedure on a cylinder taken from this tutorial:

cylinder_orig

cylinder_mod

Check out the capabilities of the package in the treatment of real CAD files in the next example:

DTMB_ffd

DTMB-5415 hull: morphing of the bulbous bow starting from an industrial .iges CAD file.

GUI

PyGeM comes with a nice GUI. For the laziest there is a video tutorial.

gui_PyGeM
PyGeM GUI: how it appears when it pops up.