Generating positioned roughness elements

Read parameters from file

from cfdmod.use_cases.roughness_gen.parameters import PositionParams
import pathlib

import pprint

pp = pprint.PrettyPrinter()

cfg_file = pathlib.Path("./fixtures/tests/roughness_gen/position_params.yaml")
cfg = PositionParams.from_file(cfg_file)

output_path = pathlib.Path("./output/roughness_gen")

pp.pprint(cfg.model_dump())

{'bounding_box': {'end': (2400.0, 1000.0, 0.0),
                  'start': (-500.0, -1000.0, 0.0)},
 'element_params': {'height': 5.0, 'width': 5.0},
 'spacing_params': {'line_offset': 10.0,
                    'offset_direction': 'x',
                    'spacing': (20.0, 20.0)},
 'surfaces': {'disk': './fixtures/tests/roughness_gen/disk/disk.lnas',
              'loft': './fixtures/tests/roughness_gen/loft/loft.lnas'}}

Read mesh files and get surfaces bounding box

from lnas import LnasFormat
import trimesh

import numpy as np

surfaces_read: dict[str, LnasFormat] = {}
surfaces_mesh: dict[str, trimesh.Trimesh] = {}

mesh_bbox = [
    (float("inf"), float("inf"), float("inf")),
    (float("-inf"), float("-inf"), float("-inf")),
]

for sfc, sfc_path in cfg.surfaces.items():
    lnas = LnasFormat.from_file(pathlib.Path(sfc_path))
    surfaces_read[sfc] = lnas
    surfaces_mesh[sfc] = trimesh.Trimesh(
        vertices=lnas.geometry.vertices, faces=lnas.geometry.triangles
    )

    min_point = surfaces_read[sfc].geometry.vertices.min(axis=0)
    max_point = surfaces_read[sfc].geometry.vertices.max(axis=0)

    if any(min_point < mesh_bbox[0]):
        mesh_bbox[0] = np.array([min_point, mesh_bbox[0]]).min(axis=0)
    if any(max_point > mesh_bbox[1]):
        mesh_bbox[1] = np.array([max_point, mesh_bbox[1]]).max(axis=0)

s = np.array([cfg.bounding_box.start, mesh_bbox[0]]).max(axis=0)
e = np.array([cfg.bounding_box.end, mesh_bbox[1]]).min(axis=0)
bbox_to_use = np.array([s, e])

bbox_to_use

array([[ -500.        , -1000.        ,   751.74133301],
       [ 2400.        ,  1000.        ,     0.        ]])

Calculate generation parameters

lx = bbox_to_use[1][0] - bbox_to_use[0][0]
ly = bbox_to_use[1][1] - bbox_to_use[0][1]

if cfg.spacing_params.offset_direction == "x":
    Nx = int((lx - cfg.spacing_params.line_offset) // (cfg.spacing_params.spacing[0]) + 1)
    Ny = int(
        (ly + cfg.spacing_params.spacing[1])
        // (cfg.element_params.width + cfg.spacing_params.spacing[1])
    )
else:
    Nx = int(lx // (cfg.spacing_params.spacing[0]) + 1)
    Ny = int(
        (ly + cfg.spacing_params.spacing[1] - cfg.spacing_params.line_offset)
        // (cfg.element_params.width + cfg.spacing_params.spacing[1])
    )

Nx, Ny, lx, ly

(145, 80, 2900.0, 2000.0)

Setup generation parameters

from cfdmod.use_cases.roughness_gen.parameters import GenerationParams

generation_params = GenerationParams(
    N_elements_x=Nx,
    N_elements_y=Ny,
    element_params=cfg.element_params,
    spacing_params=cfg.spacing_params,
)

Build elements and apply linear patterns

from cfdmod.use_cases.roughness_gen import build_single_element, linear_pattern

triangles, normals = build_single_element(generation_params.element_params)

single_line_triangles, single_line_normals = linear_pattern(
    triangles,
    normals,
    direction=generation_params.spacing_params.offset_direction,
    n_repeats=generation_params.single_line_elements,
    spacing_value=generation_params.single_line_spacing,
)

full_triangles, full_normals = linear_pattern(
    single_line_triangles,
    single_line_normals,
    direction=generation_params.perpendicular_direction,
    n_repeats=generation_params.multi_line_elements,
    spacing_value=generation_params.multi_line_spacing,
    offset_value=generation_params.spacing_params.line_offset,
)

# Offset to match bounding box limits
full_triangles[:, :, 0] += bbox_to_use[0][0]
full_triangles[:, :, 1] += bbox_to_use[0][1]

print(
    "Single element:",
    f"Vertices count: {triangles.shape[0] * triangles.shape[1]}",
    f"Triangles count: {len(triangles)}",
)
print(
    "Single line",
    f"Vertices count: {single_line_triangles.shape[0] * single_line_triangles.shape[1]}",
    f"Triangles count: {len(single_line_triangles)}",
)
print(
    "Replicated lines",
    f"Vertices count: {full_triangles.shape[0] * full_triangles.shape[1]}",
    f"Triangles count: {len(full_triangles)}",
)

Single element: Vertices count: 6 Triangles count: 2
Single line Vertices count: 870 Triangles count: 290
Replicated lines Vertices count: 69600 Triangles count: 23200

Create profiles with x normal planes

x_vals = np.unique(full_triangles[:, :, 0].reshape(-1))

profiles: dict[float, np.ndarray] = {}

for sfc_label, mesh in surfaces_mesh.items():
    x_vals_to_use = x_vals[
        (x_vals >= mesh.vertices[:, 0].min()) & (x_vals <= mesh.vertices[:, 0].max())
    ]

    for x_val in x_vals_to_use:
        section_slice = mesh.section(
            plane_origin=[x_val, 0, 0],
            plane_normal=[-1, 0, 0],
        )
        vertices = np.array(section_slice.to_dict()["vertices"])
        if x_val in profiles.keys():
            profiles[x_val] = np.concatenate((profiles[x_val], vertices))
        else:
            profiles[x_val] = vertices

for k, profile in profiles.items():
    profile = np.array(sorted(profile, key=lambda x: x[1]))

Offset geometry with surface information

for tri in full_triangles:
    left_bound = tri[:, 1].min()
    right_bound = tri[:, 1].max()
    prof = profiles[tri[0][0]][
        (profiles[tri[0][0]][:, 1] >= left_bound) & (profiles[tri[0][0]][:, 1] <= right_bound)
    ]
    while len(prof) == 0:
        value_offset = generation_params.element_params.width
        left_bound -= value_offset
        right_bound += value_offset
        prof = profiles[tri[0][0]][
            (profiles[tri[0][0]][:, 1] >= left_bound) & (profiles[tri[0][0]][:, 1] <= right_bound)
        ]

    offset_val = prof[:, 2].min()
    tri[:, 2] += offset_val

Export generated geometry

from cfdmod.api.geometry.STL import export_stl

export_stl(output_path / "positioned_elements.stl", full_triangles, full_normals)