# 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 !

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:

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

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:

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

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

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

# 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)

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

# cycle on the faces to update the control points position
# init some quantities
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)
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()
n_faces += 1
faces_explorer.Next()

writer.Write(self.outfile)```

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

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

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.

PyGeM GUI: how it appears when it pops up.