Three-dimensional ray tracing in python

I was looking for an algorithm for ray tracing in 3-D, and it took me some while to find one which can easily be programmed. Finally, I ran into this page, where a 3-D ray-tracing algorithm is clearly explained along with an example for programming. I decided to share my Python 2.7 implementation so that searchers either find this code and use it, or spend less time finding the reference to program their own code.

I wrote the algorithm in Python, slightly modified compared to the one in the link such that it is more suitable for checking whether a point is inside a (not necessarily convex) polygon.

import numpy as np

# source:

def ray_intersect_triangle(p0, p1, triangle):
    # Tests if a ray starting at point p0, in the direction
    # p1 - p0, will intersect with the triangle.
    # arguments:
    # p0, p1: numpy.ndarray, both with shape (3,) for x, y, z.
    # triangle: numpy.ndarray, shaped (3,3), with each row
    #     representing a vertex and three columns for x, y, z.
    # returns: 
    #    0.0 if ray does not intersect triangle, 
    #    1.0 if it will intersect the triangle,
    #    2.0 if starting point lies in the triangle.

    v0, v1, v2 = triangle
    u = v1 - v0
    v = v2 - v0
    normal = np.cross(u, v)

    b = np.inner(normal, p1 - p0)
    a = np.inner(normal, v0 - p0)
    # Here is the main difference with the code in the link.
    # Instead of returning if the ray is in the plane of the 
    # triangle, we set rI, the parameter at which the ray 
    # intersects the plane of the triangle, to zero so that 
    # we can later check if the starting point of the ray
    # lies on the triangle. This is important for checking 
    # if a point is inside a polygon or not.
    if (b == 0.0):
        # ray is parallel to the plane
        if a != 0.0:
            # ray is outside but parallel to the plane
            return 0
            # ray is parallel and lies in the plane
            rI = 0.0
        rI = a / b

    if rI < 0.0:
        return 0

    w = p0 + rI * (p1 - p0) - v0

    denom = np.inner(u, v) * np.inner(u, v) - \
        np.inner(u, u) * np.inner(v, v)

    si = (np.inner(u, v) * np.inner(w, v) - \
        np.inner(v, v) * np.inner(w, u)) / denom
    if (si < 0.0) | (si > 1.0):
        return 0

    ti = (np.inner(u, v) * np.inner(w, u) - \
        np.inner(u, u) * np.inner(w, v)) / denom
    if (ti < 0.0) | (si + ti > 1.0):
        return 0

    if (rI == 0.0):
        # point 0 lies ON the triangle. If checking for 
        # point inside polygon, return 2 so that the loop 
        # over triangles can stop, because it is on the 
        # polygon, thus inside.
        return 2

    return 1

I hope it is useful for those who found this page!

Curvature-based refinement in snappyHexMesh

On appealing feature of snappyHexMesh is the auto-refinement option depending on surface curvature. It is possible to define a minimal and a maximal refinement level per surface or region, and SHM will decide based on the set ‘resolveFeatureAngle’ whether to further refine that surface (up to the maximum refinement level, of course).

In my experience, however, this feature can create some grid problems. Especially at the interface between multiple refinement levels, snappyHexMesh seems to have troubles snapping the mesh correctly to the ship surface. What happens, is that if the surface locally has higher curvature, snappyHexMesh will refine the cells there. This may be limited to a single cell. Depending on the ‘nCellsBetweenLevels’ setting, more cells will be refined. However, SHM still generates an ‘island’ with an interface between larger and smaller cells around the perimeter. The initially snapped mesh may still look fine, but if the cells at the border between refinement levels have distorted shapes, adding layers is becoming more and more of a problem.

Example ship

To sketch an example, I’m trying to get a good mesh around an inland ship hull. For speed of meshing, I am currently focusing on the most complex part, the stern (an example is given in the next figure). An inland ship stern is far from smooth, due to propeller tunnels. The skirts of these tunnels have to prevent the propeller from ventilating (sucking in air from the water surface) in case the ship has to operate at low draft. One may imagine that fitting a mesh around this is not easy, but with some practice it will work out.

Stern of an inland ship. Note the tunnel skirts making meshing pretty complex.

Fitting the grid around the propeller tunnel: locally refined zones

A propeller tunnel like this is a nice candidate for using the auto-refinement function. Namely, upstream the deform smoothly into the parallel midbody of the ship, hence reducing curvature and requiring less small grid cells. The result, however, is not exactly nice:

A close-up of the mesh in the region where the tunnel blends into the parallel midbody of the ship

There are multiple small islands at places where the curvature of the surface exceeds a certain threshold locally. Again, at interfaces between different refinement levels, snappyHexMesh has problems – in my experience – fitting a nicely snapped grid to the STL surface.

Snapping artifacts at interfaces between refinement levels

In the figure above, such problems are not immediately obvious, but the following two pictures may clearify the problem.

Close-up of the snapped surface of the tunnel skirt, from the inside of the ship. Snapping problems are clearly visible.

The first picture shows the surface snapped to the STL file by snappyHexMesh. Multiple badly snapped cells can be observed. Such artifacts prevent layers from being grown at the surface.The next picture shows the mesh at the tunnel skirt. Comparing both pictures, we can see that artifacts in the snapped surface seem to occur at the interface between multiple refinement levels.

Close-up of the snapped mesh at the tunnel skirt.

My solution to get a better snapped mesh

My solution to this was to not allow snappyHexMesh to automatically refine the surface based on curvature. This means that you have to sub-divide your STL in multiple parts (the original surface on which the STL was based is in Rhino, so that was not a big problem in my case) and have the minimum and maximum refinement levels equal per surface (between surfaces, they can vary). When sub-dividing the STL, it is a good idea to ensure that the edge between separate patches is in an area of limited curvature; otherwise SHM will still have problems in those areas!

A small relieve, however, is that the problem mostly seems to occur if the interface between multiple refinement levels does not resemble a straight or smooth line. At the parallel midbody of the ship used in this post, for example, the distribution of curvature is constant: flat at the bottom, curved with radius R at the bilge, and flat at the side again (and that for 70 meters in length). In such a case, the interface between refinement levels will end up to be a straight line. SnappyHexMesh handles this without problems. But specifically for areas with a lot of curvature changes, using fixed refinement levels per surface patch seems a good idea.

A mesh check before layer addition

A nice way to check this if layers can be fitted after snapping is through the snappyHexMesh log: at the starting point of a feature snapping iteration (Morph iteration xxx), snappyHexMesh prints the minimal and maximal distance from the mesh to the STL surface. A a rule of thumb, I assume maximal distance here should be (considerably) less than the thickness of the first layer.

This applies when ‘relativeSizes’ is false. If you’re using relative layer sizing (based on the local grid cell size), note that the maximum snapping distance is a global value and if layers locally are larger due to relative layer sizing, layers may still be fitted nicely. I am using absolute layer sizing since I do not want layers to suddenly become smaller if my local patch has been refined to smaller cells (namely, the boundary layer thickness does not suddenly become different if the underlying patch is finer or coarser).

Inland ship stern optimization in shallow water


Erik Rotteveel, Auke van der Ploeg, Robert Hekkenberg


June 21, 2017 in Ocean Engineering: Special issue on Marine Applications in Fluid Dynamics


Inland ships continuously operate in restricted waters, where the depth and width are regularly less than twice the ship’s draft and four times ship breadth, respectively. In restricted water, the flow around the hull changes compared to that in unrestricted water due to presence of the fairway bottom and sides, that lead to increased return flow, stronger squat effects and changes in the wave pattern produced by the ship. If these changes to the flow are significant, it is worthwhile to optimize the hull form for shallow or confined water rather than for unrestricted water. This paper specifically focuses on the effects of water depth on inland ship stern optimization. It presents the optimization of propulsion power for various water depths using a parametric inland ship stern shape, CFD and surrogate modeling. The change of parameter influence in different water depths is analyzed and explained by means of flow visualization. Using Pareto fronts, a trade-off is shown: propulsion power in shallow water can be decreased at the cost of increased propulsion power in deep water and vice versa.

Full paper

Please visit to read the full paper.

A ducted actuator disk

This is my first post on this site, since I only recently started practicing with OpenFOAM. So far, however, I noticed that for an open-source software suite that is not thát widely used, it is quite mature. However, not all settings are well documented, and unfortunately a lot of users of the cfd-online forums don’t post their solutions if they got their case finally working. That is also why I started this whole website, to show you how I approached certain cases that are relevant to my own work.

Download the case files here: duct_case.tar

A ducted propeller in OpenFOAM

In my PhD project, I focused on inland ships. These are self-propelled ships and due to limitations on their dimensions (locks, bridges, fairway dimensions), these ships are really full (that is, a high block coefficient C_B = \nabla / (L \cdot B \cdot T)). Together with limited draft imposed by limited water depth, these ships can only fit a small propeller. To increase the amount of thrust and power per square meter of propeller area, and to improve efficiency, inland ships are equipped with ducted propellers.

A 19A duct

The fixedJump boundary condition

Looking around the cfd-online forums and OpenFOAM-related websites in general, it seemed like modeling an actuator disk as used in Maritime Hydrodynamics, a pressure jump, is not straightforward. It apparently is not possible to – without coding – add body forces in a certain region of the mesh (at least I did not manage to, if you know how, leave a comment!). Looking through the C++ guide of OpenFOAM 4.1, I came across the FixedJump boundary condition. It uses a cyclic boundary condition at its base, but allows a fixed difference for a certain field variable between the two connected boundaries. This is exactly what an actuator disk requires. All field variables should be continuous, except for a fixed jump in pressure. The value of the jump can be computed from the required total thrust force and the area of the disk.

Creating baffles in snappyHexMesh

In order to implement the fixedJump boundary condition, we require that there are two patches, with the same points and topology and preferably very, very close to each other. I tried multiple approaches – I did not find the right solution on the cfd-online forum, such as creating a closed solid inside the duct, and using the opposing sides of this solid to implement the cyclic and fixedJump boundary conditions. However, snappyHexMesh puts different, non-conformal vertices on both sides. The cyclicAMI boundary conditions – which should allow for a mapping between two patches with unequal topolgy, did not work either.

Then I found: OpenFOAM 2.2.0: snappyHexMesh. Within the refinementSurfaces sub-dictionary in snappyHexMeshDict, it is possible to define a faceZone. This faceZone can be defined as a baffle, i.e. a very thin surface. With this setting, snappyHexMesh will ensure that interior grid points are snapped to this surface, and finally create two patches (a master and a slave patch). In my case, the surface got the name ‘disk1’, and therefore two cyclic patches, ‘disk1’ and ‘disk1_slave’, will turn up in the resulting mesh.

My approach

I used the creation of baffles in snappyHexMesh as well as the fixedJump boundary condition to create the pressure jump across the ducted propeller. Furthermore, I used simpleFoam to solve the flow equations. simpleFoam is fast and thus suitable to find out how a ducted propeller should be modeled.

Adding the faceZone

First, I have to add the faceZone to the snappyHexMeshDict file. This is done as shown in the following code block. During meshing, snappyHexMesh will create two internal patches: disk and disk_slave. The reason it uses ‘disk’ is because it is defined with that name in the geometry section of the snappyHexMeshDict file.

Unfortunately, I did not manage to find how to directly set the internal patches as cyclic patches from snappyHexMesh, it seems to ignore additional patchInfo settings for faceZones. However, this can be overcome using createPatch (or createBaffles, if you like).


        { file "disk1.eMesh"; level 1; }

        // add the duct (and other surfaces) here
            level (1 1);
            faceZone diskFaces;
            faceType baffle;


    // Remaining settings

Some snapping and layer addition experience

In the case that can be downloaded at the bottom of the page, note that I use two separate dictionaries for snappyHexMesh: one for refinement and snapping, and a second one to do the layer addition. This allows me to first check whether the surface is snapped smoothly to the duct (which is fortunately is), before adding layers (which otherwise fails in places where the mesh is badly snapped and ‘bumps’ occur on the patches).

Apart from the separated snappyHexMeshDict files, I also subdivided the duct STL in four parts. This allows me to have the leading and trailing edge separated, and simply fix the refinement level of the mesh on both. The reason for this is that if snappyHexMesh is refining based on feature angles or local curvature, small ‘islands’ of higher refinement levels occur, and those seem to be locations where snappyHexMesh clearly has troubles getting a well-snapped mesh. Hence, I use equal values for minimum and maximum refinement per surface. The border between the lower and  higher refinement levels then becomes an almost straight line, something that snappyHexMesh seems to work well with.

The snapped mesh around the duct. Note the (almost) sharp border between different refinement levels. This seems to be important for snappyHexMesh to snap the surface well enough. I did need to have 2 buffer layers (nCellsBetweenLevels) though. With this setting at 1, there would still be some bumps.

Another thing I found out for this specific case is that in the snapControls sub-dictionary, it is useful to keep nSmoothPatch low, but increasing nSolveIter doesn’t hurt and may increase the quality of the final mesh. Futhermore, although I added a specific feature edge for the connection between the duct and the disk patch (that will later be used to put the pressure jump on), the setting explicitFeatureSnap does not make things better. In fact, it was better for me to use implicitFeatureSnap.

When adding layers, I feel that snappyHexMesh prefers to have each of the surfaces on which you want layers specified separately rather than just using the patch group ‘duct’ (which I in fact use to define the no-slip BC later on). Hence, like this:

  leadingedge	{ nSurfaceLayers 3; }
  trailingedge	{ nSurfaceLayers 3; }
  innerring	{ nSurfaceLayers 3; }
  outerring	{ nSurfaceLayers 3; }

With regard to settings, increasing nSmoothNormals – which smoothens the normals of the interior grid, not those on the patch – works pretty well. Also in cases with a more complex geometry, I found that it is worthwhile to wait a bit for this process. With nSmoothSurfaceNormals, however, you should be careful. For example, the normals close to the trailing edge (but still on the ‘flat’ part) will start pointing backward if too many surface smoothing iterations are used, which creates a distorted layer mesh.

Adding the fixedJump boundary condition

With a nicely snapped mesh including good-looking layers (y+ is 100, which is OK when using wall functions), the boundary conditions on the patches can be defined. In the case you can download, you’ll see that createPatch is used to ensure that all patch names are correct, that all parts of the duct are grouped under ‘duct’, and – most important – that the internal baffle faces are set to be cyclic patches. I also give them somewhat better names, namely disk_in and disk_out.

Next, the boundary conditions need to be defined. This is done in the ‘0’ folder. We need to implement a pressure jump, so I’ll start with what is required for the ‘p’ file. I do assume, however, that other boundary conditions for ‘p’ are known to  you already.

    type	fixedJump;
    patchType	cyclic;
    jump	uniform -0.8;
    type	fixedJump;
    patchType	cyclic;

Note that in the above code, only one of the patches gets assigned the jump! This is actually the master patch. Remember that snappyHexMesh created both ‘disk’ and ‘disk_slave’. Using createPatch, I renamed the first to ‘disk_in’ and the latter to ‘disk_out’. Hence, the jump should be defined on ‘disk_in’, which is the master patch.

As for the value of the jump, I simply wanted a thrust coefficient (of the propeller only!) between 1.0 and 2.0:

C_T=\frac{T}{0.5 \cdot \rho \cdot V^2 \cdot A_0}

I used a thrust coefficient of 1.6, and since the used duct has a radius of 1.0 meter, and I use a speed of 1.0 meter as well, the total thrust is approximately 2500 newtons. Dividing this by the are of the disk (which conveniently is exactly \pi), the pressure jump is around 800 newtons per square meter. However, simpleFoam, being an incompressible, single-phase solver, uses the pressure divided by density for the ‘p’ variable! Thus, to have a 800 pascal pressure jump, the jump in the above code must be set to 0.8. Finally, it is negative because of the orientation I used: the flow is directed from large to small x.

For other variables than pressure, the boundary condition on the disk patches is simply cyclic, as displayed in the following:

    type cyclic; 
    type cyclic;

After setting the boundary conditions, I used a k\epsilon turbulence model with wall functions, so that a very fine grid around the duct is not required. Further settings are pretty standard and are taken from the tutorials that are shipped with openFoam.

Note, however, that the grid downstream of the duct should be finer in order to really capture the flow phenomena occuring there. At this moment, the (circular) wake of the duct is quickly smoothened out across the coarse grid cells downstream of the duct. However, this is just a test case and me practicing with OpenFOAM.

Download the case files

I hope this post helped for those needing a ducted actuator disk in OpenFOAM. The files can be downloaded through this link: duct_case.tar