Yours truly is an avid gamer, @fadoobaba (YouTube). The ray casting algorithm๐ธ; a fundamental ๐งฑ technique used in video game ๐ฎ development and computer graphics, has been implemented within the finite difference method ๐ code yours truly has been developing. In abundant spare time, of course๐ผ. In this post, this method is explained. For details about ray casting using matplotlib.path, refer to here. For validation of the code, refer to here, backwards-facing step, curved boundaries, here (moving cylinder) and here.
NOTE: This method requires a GPU, if dear readers don't have a GPU then please stop being peasants... ๐
If you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2025). S-IBM (https://fluiddynamicscomputer.blogspot.com/2025/07/a-simplified-immersed-boundary-method.html), Blogger. Retrieved Month Date, Year
The first step is to setup the polygon ๐ . The polygon = cp.column_stack((x1, y1)) statement is used to combine x1 and y1 into a single array representing the polygon vertices. The following statements are used to store the x and y coordinates of the polygon vertices and nvert is the total number of vertices.
px = polygon[:, 0]
py = polygon[:, 1]
nvert = len(polygon)
The boolean masks ๐บ are then initialized. The following arrays track whether a grid point is inside the polygon based on horizontal and vertical ray intersections. False = outside ❌, True = inside ✅. Then a loop is implemented to for each edge ๐ of the polygon, defined by vertices i (current) and j (previous), with % nvert ensuring the polygon is closed (last vertex connects back to the first). (xi, yi) and (xj, yj) define the current edge and previous edge.
A check ⁉️ is performed to ensure a ray crosses the polygon edge once. cond1 statement checks if the test point lies between the y-values of the edge's endpoints i.e., a ray could cross it. intersect_x finds where the edge crosses a horizontal line at test_y. cond2 checks if the intersection lies to the right of the test point. ^= is XOR toggles the "inside" state each time the ray crosses an edge. The addition of small term ~1e-16 prevents division by zero for horizontal edges. Similar method is applied to verify the points using vertical ray. cond3 checks if the edge crosses a vertical ray (top to bottom) from (test_x, test_y). slope1 and intersect_y computes where the vertical ray at x=test_x intersects the edge. A point is inside only if both horizontal and vertical rays classify it as inside. A point is considered inside the foil only if both ray checks are true. curve is the 2D boolean mask of grid points inside the foil body.
cond1 = ((yi > test_y) != (yj > test_y))
slope = (xj - xi) / (yj - yi + 1e-16)
intersect_x = slope * (test_y - yi) + xi
cond2 = test_x < intersect_x
horizontal_inside ^= cond1 & cond2
cond3 = ((xi > test_x) != (xj > test_x))
slope1 = (yj - yi) / (xj - xi + 1e-16)
intersect_y = slope1 * (test_x - xi) + yi
cond4 = test_y < intersect_y
vertical_inside ^= cond3 & cond4
inside = horizontal_inside & vertical_inside
curve = inside.reshape(X.T.shape)
Using both horizontal and vertical rays reduces false positives (e.g., near sharp corners). The & ensures only points unambiguously inside are marked. Within Fig. 1, the points on the shape boundary, points inside ๐ชฐ and outside ๐ the shape boundary are shown. A point is a boundary if it is inside the body, but any of its neighbors is outside. Following statements are used to mark the boundary of the object. Boundary detects the "skin" of the body for applying no-slip, force, or stress conditions.
interior = curve[1:-1, 1:-1]
right = curve[2:, 1:-1]
left = curve[:-2, 1:-1]
top = curve[1:-1, 2:]
bottom = curve[1:-1, :-2]
boundary = interior & ~(right & left & top & bottom)
The statements boundary_indices = cp.where(boundary_mask) ... valid_right = (boundary_indices[0] + 1 < X.shape[0]) & (~curve[right_neighbors]) ... etc. extract boundary indices and their valid neighbors i.e. these lines extract the (i, j) grid indices of the surface and locate which neighbor cells are valid fluid neighbors (outside the body and within domain bounds). These are needed to compute normals or apply boundary conditions via interpolation or extrapolation from the fluid.
Fig. 1, Mesh cells
For validation of the results from present simulations, the case of flow around a circular cylinder is selected. Fig. 1 shows the results from the code at Re 200. The drag coefficient obtained from this code is 1.396 while from the literature, the value is at 1.4 [1]. The lift coefficient from the code is at 0.000134. Within Fig. 2, top row shows u and v velocity components and bottom row shows pressure and vorticity.
Fig. 2, The flow-field
For the second an more complex validation case, a swimming fish is simulated. The lift and drag coefficients obtained from the simulation are compared with experimental results. The drag coefficient from the current code is at 0.342 while from the published literature, the value is at 0.348. Maximum lift coefficient is at 7.243 and 8 from the present code as compared to the published literature [2]. Within Fig. 3 the u and v velocities, pressure and velocity streamlines are shown for St = 0.8 and Reynolds number of 500. The computational mesh near the fish is shown in Fig. 4. Within Fig.4, a zoomed in view towards the right shows the mesh at the trailing edge of the fish.
Fig.3, Post-processing of results
Fig. 4, The mesh
This method allows handling arbitrary deforming / non-deforming shapes on a fixed Cartesian grid. In summary, the method has the following steps.
3. Use ray casting to check if points are inside shape
4. Build a boolean mask of body region
5. Identify surface (boundary) points
6. Extract boundary indices for physics coupling
References
[1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014
[2] Fulong Shi, Jianjian Xin, Chuanzhong Ou, Zhiwei Li, Xing Chang, Ling Wan; Effects of the Reynolds number and attack angle on wake dynamics of fish swimming in oblique flows. Physics of Fluids 1 February 2025; 37 (2): 025205 doi.org/10.1063/5.0252506
Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!
In this post, details will be shared about a simplified version of the Immersed ๐คฟ Boundary Method for fluid dynamics ๐จ simulations which yours truly has developed in abundant spare time ⏲, referred to as S-IBM ๐ค. The purpose of this code is for teaching ๐ purposes only. The mesh is Cartesian and the method used has no measures in place to resolve the stairstep ๐ mesh. The syntax share in this post is for Python ๐, but this code can be translated to other languages such as C or MATLAB, easily.
The process starts with creating a polygon, a circular cylinder ⭕ for this post. A polygon is created using parametric equations for the circle. NOTE: Before writing these lines, several parameters such as domain size, mesh size, number of points in each coordinate direction etc. have been already defined and calculated. For more information, read here.
theta = np.linspace(0, 2 * np.pi, 100) # create equally spaced angles from 0 to 2 * pi
Path from matplotlib is used to create a polygon approximation of a circle or an airfoil by using the x and y co-ordinates of the circle (x1, y1) or any shape. Next, coordinate pairs of the grid points are created using vstack and ravel commands. A check is performed to verify if a grid point falls inside the polygon. custom_curve is a boolean array. The values inside this array are True only for the points inside the polygon. Transpose is used to keep the curve orientation consistent with the Navier-Stokes equations.
Another boolean mask is created for boundary points using the boundary_mask statement. Points that are inside the polygon (circle) are marked using using the interior statement.
boundary_mask = np.zeros_like(curve, dtype=bool)
interior = curve[1:-1, 1:-1]
Next a check is performed to mark each interior point's immediate neighbor. right means point to the right of each interior point, bottom means point to the bottom of each interior point etc. boundary statement is used to identify the boundary points (fluid-solid interface) according to the following logic. The point itself is inside the circle (interior = True) and at least one neighbor is outside. boundary_mask[1:-1, 1:-1] is used to fill the boolean with boundary. Consider a 5×5 a grid (zoomed near the polygon boundary). X marks solid points (inside circle), O marks fluid points, as shown in Fig. 1. For example, if interior = True, right = True, left = True, top = True and bottom = True → boundary = True & ~(True) = False i.e. an interior point. And for example, interior = True, right = False, left = True, top = True, bottom = True → boundary = True & ~(False) = True, i.e. a point on the right boundary. For points outside the polygon, interior = False → boundary = False (regardless of neighbors).
Fig. 1, A grid example
right = curve[2:, 1:-1]
left = curve[:-2, 1:-1]
top = curve[1:-1, 2:]
bottom = curve[1:-1, :-2]
boundary = interior & ~(right & left & top & bottom)
boundary_mask[1:-1, 1:-1] = boundary
custom_curve_boundary = boundary_mask
In the next step, boundary_indices obtains indices of all boundary points. right_neighbors, left_neighbors etc. calculates indices of neighboring points in all directions. valid_right, valid_left etc. checks if the neighbor is within domain bounds (not at edge) and is outside the circle (fluid point). This allows proper application of no-slip boundary conditions. If for example, valid_right is True, pressure is taken from the right neighbor to be applied to the corresponding boundary index. A visualization of the points is shown in Fig. 2.
Fig. 2, Depiction of grid points in various regimes
While plotting Fig. 2, several points are skipped. This is done to make plot clearer. This code can be easily implemented in other languages such as MATLAB (inpolygon), C (ray-casting) etc. At Re 200, this code predicts a Drag coefficient of 1.395 and a Lift coefficient of 0.053. These values are well within rage of the values reported by [1]. In this code, there is no interpolation of scheme or smoothing employed. The lift and drag force coefficients for one shedding cycle are shown in Fig. 3. Within Fig. 4, u and v velocities (top left and right), pressure and streamlines (bottom left and right) are shown.
Fig. 3, red = Cd and black = Cl
Fig. 4, Post processing from the simulations.
References
[1] Braza M, Chassaing P, Minh HH. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder. Journal of Fluid Mechanics. 1986;165:79-130. doi:10.1017/S0022112086003014
Thank you for reading! If you want to hire me as a post-doc researcher in the fields of thermo-fluids and / or fracture mechanics, do reach out!
This is probably the final step in the series. Unless, one fine morning yours truly randomly decides to write about moving and morphing airfoils ๐ i.e. a true virtual wind tunnel ๐ค.
This whole series is a continuation of a series created by Dr. Barbra ๐ฉ๐ซ to teach about coding ๐ณ of the Navier-Stokes ๐ฌ equations. The code made available this blog is an independent continuation of her open-source work. Yours truly has coded the scaled version of the equations with a modified pressure Poisson equation. The original series ends with the example of the channel flow ๐. So far, yours truly has added flow around obstacles both curved and non-curved; flow around moving obstacles and species transport to the series. The Fluid-Structure Interaction (FSI) version of the code ๐ป is available here.
This post is about flow around shape changing ๐ชฑ obstacles i.e. morphing wings. The development of these codes shared here, including the one mentioned in this post, started in the year 2020 and continued till the fall of 2023. All this effort was to create a digital CV so that yours truly can get himself accepted in to a funded PhD program. That happened in fall 2023 ๐.
The frequency of oscillations, the equations for amplitude and lateral oscillations 〰️ are all taken from the published literature about swimming fish ๐. As a starting point, the NACA 0012 airfoil is used. A polygon is then created. Each grid index is then tested for presence inside or outside the polygon. An index is marked as the boundary of the polygon if one index has anyone neighboring index outside the polygon. The neighboring boundary indices are used to apply Neumann boundary conditions for pressure. The boundary indices are used to apply no-slip Dirichlet boundary conditions on the polygon boundary. More details can be read here.
Because the camber of airfoil changes at each time-step, approximately 400 points are required on the airfoil boundary to properly capture the time-varying camber of the airfoil i.e. the fish ๐. The output from the code is shown in Fig. 1. Within Fig. 1, top row has u and v components of velocity. Meanwhile, the bottom row shows the pressure and streamlines. The animation is shown in Fig. 2.
If you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2025). 18th Step of the 12 Steps to Navier-Stokes (https://fluiddynamicscomputer.blogspot.com/2025/05/18th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year
Fig. 1, The result from code
Fig. 2, Animation
Code
#Copyright <2025> <FAHAD BUTT>
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
y1 = foil_center_y + np.concatenate([y_top, y_bottom[::-1]]) # foil y
foil_path = Path(np.column_stack((x1, y1))) # foil region (polygon)
points = np.vstack((X.ravel(), Y.ravel())).T # check if inside region
custom_curve = foil_path.contains_points(points).reshape(X.shape) # test all grid points for inclusion in airfoil, contains_points returns True for each point inside the airfoil (polygon)
curve = custom_curve.T # transpose to match grid orientation
right = curve[2:, 1:-1] # boolean arrays indicating if neighbors of each interior point are inside airfoil (polygon)
left = curve[:-2, 1:-1]
top = curve[1:-1, 2:]
bottom = curve[1:-1, :-2]
boundary = interior & ~(right & left & top & bottom) # mark if point is inside foil and at least one neighbor is outside foil (mark interface between solid and fluid)
boundary_mask[1:-1, 1:-1] = boundary # fill the empty boundary
custom_curve_boundary = boundary_mask.T # transpose back
boundary_indices = np.where(custom_curve_boundary.T) # mark boundary indices
right_neighbors = (boundary_indices[0] + 1, boundary_indices[1]) # right index (neighboring grid points in each direction)
left_neighbors = (boundary_indices[0] - 1, boundary_indices[1]) # left index
top_neighbors = (boundary_indices[0], boundary_indices[1] + 1) # top index
bottom_neighbors = (boundary_indices[0], boundary_indices[1] - 1) # bottom index
# identify valid neighbor indices, ensure the neighbor is not inside airfoil (polygon) i.e. is inside fluid)
FOANSS, i.e. Fadoobaba's Open Advanced Navier-Stokes ๐ฌ Solver is now capable to simulate flow with moving objects ๐ค. In this example, the objects are square and circular cylinders. These blog posts could have been a master's thesis or a peer-reviewed publication but yours truly decided to are write here for educational purposes! The
The simulated case is called the case of flow around an obstacle. For validation, read. This code and blog post is not endorsed or approved by Dr. Barba, I just continue the open-source work of her. The 13th, 14th, 15th and 16th Steps are available for reading along with code. The resulting motion of the square cylinder is shown in Fig. 1. The Fig. 1 is colored using the horizontal velocity. The Fig. 2 shows the example of moving circle. Within Fig. 2, the results of u and v velocities are shown on top while the pressure and streamlines are shown at the bottom. The code for stationery curved / arbitrary boundaries is made available here. A customary remainder that you are reading a blog, not a Q1 journal ๐.
If you plan to use these codes in your scholarly work, do cite this blog as:
Fahad Butt (2025). 17th Step of the 12 Steps to Navier-Stokes (https://fluiddynamicscomputer.blogspot.com/2025/04/17th-step-of-12-steps-to-navier-stokes.html), Blogger. Retrieved Month Date, Year
Fig. 1, The animation
Code
#%%
#Copyright <2025> <FAHAD BUTT>
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
plt.contourf(X, Y, u1.T, 128, cmap = 'jet') # plot contours
# plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
plt.figure(dpi = 200)
plt.contourf(X, Y, v1.T, 128, cmap = 'jet') # plot contours
# plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
plt.figure(dpi = 200)
plt.contourf(X, Y, p1.T, 128, cmap = 'jet') # plot contours
# plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
plt.figure(dpi = 200)
plt.streamplot(X, Y, u1.T, v1.T, color = 'black', cmap = 'jet', density = 2, linewidth = 0.5, arrowstyle='->', arrowsize = 1) # plot streamlines
plt.gca().set_aspect('equal', adjustable='box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.show()
Code (Curved / Slanted Boundaries)
#Copyright <2025> <FAHAD BUTT>
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
