Help Converting For Loop in Python to Julia (Indexing in For Loop)

New to Julia
1

Hi All,

I am trying to convert some relatively slow Python code to Julia. The code calculates the skewness of mesh elements for mixed meshes (quads and tris) based on the nodal coordinates. An example of what the code looks like in Python is below (element_coords is a 3D array that defines the nodal coordinates of each element):

import numpy as np

element_coords = np.zeros((2, 3, 2))
element_coords[0,:,:] = np.array([[0, 0], [1, 1], [0.5, 2]])
element_coords[1,:,:] = np.array([[2, 2], [4, 5], [3, 6]])

element_equiangular = np.ones(len(element_coords[:,0,0])) * 60
element_skewness = np.zeros_like(element_equiangular)

for jj, coords in enumerate(element_coords):
    maxAng = 0.
    minAng = (np.size(coords, 0) - 2)*180
    for ii in range(len(coords)):
        # Defining the vector from the start point to the reference point
        x1 = coords[ii, 0] - coords[ii-1, 0]
        y1 = coords[ii, 1] - coords[ii-1, 1]
        # Defining the vector from the point before the reference point to the reference point
        x2 = coords[ii-2, 0] - coords[ii-1, 0]
        y2 = coords[ii-2, 1] - coords[ii-1, 1]
        # Calculating the angle between the vectors
        ang = np.degrees(np.arctan2((x1*y2 - y1*x2), (x1*x2 + y1*y2))) #arctan2(cross(vec1, vec2), dot(vec1, vec2))
        if ang >= maxAng:
            maxAng = ang
        elif ang <= minAng:
            minAng = ang
    element_skewness[jj] = max(((maxAng - element_equiangular[jj])/(180-element_equiangular[jj])), ((element_equiangular[jj]-minAng)/element_equiangular[jj]))

I’ve tried converting this code to Julia but I am stuck on the how to convert the inner for loop to Julia. For Python when ii = 0 (first iteration), ii-1 and ii-2 represent the last and second last element of the array respectively but in Julia at the first iteration you can’t index into the array using 0 and -1 (that I know of). Any tips on how to convert this inner loop to Julia would be much appreciated, below is the start of converting this Python code to Julia:

element_coords = zeros(3, 2, 2)
element_coords[:,:,1] = [0. 0.
                         1. 1.
                         0.5 2.]
element_coords[:,:,2] = [2. 2.
                         4. 5.
                         3. 6.]
element_equiangular = ones(size(element_coords, 3)) .* 60
element_skewness = zero(element_equiangular)

function getSkew!(ele_coords, ele_equi, ele_skew)
    for jj = 1:size(ele_coords, 3)
        coords = ele_coords[:, :, jj]
        maxAng = 0.
        minAng = (size(coords,1) - 2) * 180
        for ii in 1:size(coords, 1)
            x1 = coords[ii, 1] - coords[ii-1, 1] # Doesn't work since 0 is not a valid index
            y1 = coords[ii, 2] - coords[ii-1, 2] # Doesn't work since 0 is not a valid index
            x2 = coords[ii-2, 1] - coords[ii-1, 1] # Doesn't work since -1 and 0 are not a valid index
            y2 = coords[ii-2, 2] - coords[ii-1, 2] # Doesn't work since -1 and 0 are not a valid index
            ang = rad2deg(atan((x1*y2 - y1*x2), (x1*x2 + y1*y2)))
            if ang >= maxAng
                maxAng = ang
            elseif ang <= minAng
                minAng = ang
            end
        end
        ele_skew[jj] = max((maxAng - ele_equi[jj])/(180-ele_equi[jj]), (ele_equi[jj] - minAng)/(ele_equi[jj]))
    end
end
2

I would save the length of the dimension in a variable n and then replace the indices with mod1(ii - 1, n), that should handle the wraparound that you want.

3

Easier to make use of end within the indexing.

julia> x = 1:5
1:5

julia> x[mod1(1, end)]
1

julia> x[mod1(0, end)]
5

julia> x[mod1(-1, end)]
4
2 Likes
4

Oh true I always forget end works anywhere inside the brackets.