Generates multiple layers of particles to use as wall boundary condition on DualSPHysics

generateMultiLayerWall.py

import numpy as np
from stl import mesh
import math
import numba

@numba.njit
def ray_triangle_intersection(ray_start, ray_vec, triangle):
    """Moeller–Trumbore intersection algorithm.

    Parameters
    ----------
    ray_start : np.ndarray
        Length three numpy array representing start of point.

    ray_vec : np.ndarray
        Direction of the ray.

    triangle : np.ndarray
        ``3 x 3`` numpy array containing the three vertices of a
        triangle.

    Returns
    -------
    bool
        ``True`` when there is an intersection.

    tuple
        Length three tuple containing the distance ``t``, and the
        intersection in unit triangle ``u``, ``v`` coordinates.  When
        there is no intersection, these values will be:
        ``[np.nan, np.nan, np.nan]``

    """
    # define a null intersection
    null_inter = np.array([np.nan, np.nan, np.nan])

    # break down triangle into the individual points
    v1, v2, v3 = triangle
    eps = 0.000001

    # compute edges
    edge1 = v2 - v1
    edge2 = v3 - v1
    pvec = np.cross(ray_vec, edge2)
    det = edge1.dot(pvec)

    if abs(det) < eps:  # no intersection
        return False, null_inter
    inv_det = 1.0 / det
    tvec = ray_start - v1
    u = tvec.dot(pvec) * inv_det

    if u < 0.0 or u > 1.0:  # if not intersection
        return False, null_inter

    qvec = np.cross(tvec, edge1)
    v = ray_vec.dot(qvec) * inv_det
    if v < 0.0 or u + v > 1.0:  # if not intersection
        return False, null_inter

    t = edge2.dot(qvec) * inv_det
    if t < eps:
        return False, null_inter

    return True, np.array([t, u, v])

def main():

    # Latice step 
    d = 0.15
    # Latice dimensions
    p_min = np.array([-36.001, -5.001, 50.001])
    p_max = np.array([70, 6, 101])
    nx, ny, nz = [int(math.ceil(l/d)) for l in list(p_max - p_min)]
    p_max[0] = p_min[0] + nx * d
    p_max[1] = p_min[1] + ny * d
    p_max[2] = p_min[2] + nz * d

    print(nx, ny, nz)
    max_depth = 3
        
    offsets = []
    for ox in range(-1, 2):
        for oy in range(-1, 2):
            for oz in range(-1, 2):
                if ox == 0 and oy == 0 and oz == 0:
                    continue
                offsets.append((ox,oy,oz))

    stl = mesh.Mesh.from_file('geo_blender.stl')
    print(len(stl.points))
    print(len(stl.normals))
    print(stl.points[0])
    print(stl.normals[0])

    all_hits = {}
    ray_vec = np.array([0, 0, 1.0])
    nodes = {}

    print("Collecting hits...")
    for ti, tpoints in enumerate(stl.points):
        if ti % 1000 == 0:
            print(f"  {ti} of {len(stl.points)}")
        # Calculate bounds
        min_x = min(tpoints[0], tpoints[3], tpoints[6])
        max_x = max(tpoints[0], tpoints[3], tpoints[6])
        min_y = min(tpoints[1], tpoints[4], tpoints[7])
        max_y = max(tpoints[1], tpoints[4], tpoints[7])        
        nx0 = int(np.floor((min_x - p_min[0]) / d))
        nx1 = int(np.ceil((max_x - p_min[0]) / d))
        ny0 = int(np.floor((min_y - p_min[1]) / d))
        ny1 = int(np.ceil((max_y - p_min[1]) / d))
        for ix in range(nx0, nx1 + 1):
            for iy in range(ny0, ny1 + 1):
                ray_start = p_min + np.array([ix * (p_max[0] - p_min[0]) / nx, iy * (p_max[1] - p_min[1]) / ny, 0])

                hit, intersection = ray_triangle_intersection(ray_start, ray_vec, np.float64(np.reshape(tpoints, (3,3))))
                if not hit:
                    continue
                hits = all_hits.get((ix,iy), [])
                hits.append((intersection[0], ti))
                all_hits[(ix, iy)] = hits

    nodes = np.full((nx+1, ny+1, nz+1), -1)
    for ix in range(nx + 1):
        print(f"Tracking in/out row {ix} of {nx+1}")
        for iy in range(ny + 1):
            hits = all_hits.get((ix,iy), [])
            if len(hits) == 0:
                continue
            hits.sort(key=lambda x: x[0])
            hitdist, ti = hits.pop(0)
            
            isIn = False
            for iz in range(nz + 1):
                l = iz * (p_max[2] - p_min[2]) / nz
                if l > hitdist:
                    isIn = not isIn
                    if len(hits) > 0:
                        hitdist, ti = hits.pop(0)
                    else:
                        hitdist = 1e6

                if isIn:
                    nodes[ix,iy,iz] = 0

    print("Floodfill 0")
    for ix in range(1, nx):    
        for iy in range(1, ny):
            for iz in range(1, nz):
                if nodes[ix,iy,iz] > -1:
                    continue
                for ox, oy, oz in offsets:
                    if nodes[ix+ox,iy+oy,iz+oz] == 0:
                        nodes[ix+ox,iy+oy,iz+oz] = 1

    for depth in range(1, max_depth):
        print(f"Floodfill {depth}")
        for ix in range(1, nx):    
            for iy in range(1, ny):
                for iz in range(1, nz):
                    if nodes[ix,iy,iz] != depth:
                        continue
                    for ox, oy, oz in offsets:
                        if nodes[ix+ox,iy+oy,iz+oz] == 0:
                            nodes[ix+ox,iy+oy,iz+oz] = depth+1

    with open("nodes.csv", "w") as f:
        f.write("x,y,z,d\n")
        for ix in range(nx + 1):
            for iy in range(ny + 1):
                for iz in range(nz + 1):
                    if 1 <= nodes[ix,iy,iz] <= max_depth:
                        d = nodes[ix,iy,iz]
                        x = p_min[0] + (p_max[0] - p_min[0]) * ix / nx
                        y = p_min[1] + (p_max[1] - p_min[1]) * iy / ny
                        z = p_min[2] + (p_max[2] - p_min[2]) * iz / nz
                        f.write(f"{x},{y},{z},{d}\n")

if __name__ == "__main__":
    main()