diff --git a/examples/structural_mechanics/crash/README.md b/examples/structural_mechanics/crash/README.md index 20fd191e3b..1475fa393f 100644 --- a/examples/structural_mechanics/crash/README.md +++ b/examples/structural_mechanics/crash/README.md @@ -121,6 +121,68 @@ This will install: - lasso-python (for LS-DYNA file parsing), - torch_geometric and torch_scatter (for GNN operations), +## Data Preprocessing + +`PhysicsNeMo` has a related project to help with data processing, called +[PhysicsNeMo-Curator](https://github.com/NVIDIA/physicsnemo-curator). +Using `PhysicsNeMo-Curator`, crash simulation data from LS-DYNA can be processed into training-ready formats easily. + +Currently, this can be used to preprocess d3plot files into VTP. + +### Quick Start + +Install PhysicsNeMo-Curator following +[these instructions](https://github.com/NVIDIA/physicsnemo-curator?tab=readme-ov-file#installation-and-usage). + +Process your LS-DYNA data: + +```bash +export PYTHONPATH=$PYTHONPATH:examples && +physicsnemo-curator-etl \ + --config-dir=examples/config \ + --config-name=crash_etl \ + etl.source.input_dir=/data/crash_sims/ \ + etl.sink.output_dir=/data/crash_processed_vtp/ \ + etl.processing.num_processes=4 +``` + +This will process all LS-DYNA runs in `/data/crash_sims/` and output VTP files to `/data/crash_processed_vtp/`. + +### Input Data Structure + +The Curator expects your LS-DYNA data organized as: + +``` +crash_sims/ +├── Run100/ +│ ├── d3plot # Required: binary mesh/displacement data +│ └── run100.k # Optional: part thickness definitions +├── Run101/ +│ ├── d3plot +│ └── run101.k +└── ... +``` + +### Output Formats + +#### VTP Format (Recommended for this example) + +Produces single VTP file per run with all timesteps as displacement fields: + +``` +crash_processed_vtp/ +├── Run100.vtp +├── Run101.vtp +└── ... +``` + +Each VTP contains: +- Reference coordinates at t=0 +- Displacement fields: `displacement_t0.000`, `displacement_t0.005`, etc. +- Node thickness values + +This format is directly compatible with the VTP reader in this example. + ## Training Training is managed via Hydra configurations located in conf/. diff --git a/examples/structural_mechanics/crash/tests/test_vtp_reader.py b/examples/structural_mechanics/crash/tests/test_vtp_reader.py new file mode 100644 index 0000000000..22b20646f9 --- /dev/null +++ b/examples/structural_mechanics/crash/tests/test_vtp_reader.py @@ -0,0 +1,222 @@ +# SPDX-FileCopyrightText: Copyright (c) 2023 - 2025 NVIDIA CORPORATION & AFFILIATES. +# SPDX-FileCopyrightText: All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import tempfile +import numpy as np +import pyvista as pv +import pytest +from pathlib import Path + +# Import functions from vtp_reader +import sys + +sys.path.insert(0, str(Path(__file__).parent.parent)) +from vtp_reader import ( + load_vtp_file, + extract_mesh_connectivity_from_polydata, + build_edges_from_mesh_connectivity, +) + + +@pytest.fixture +def simple_vtp_file(): + """Create a simple VTP file for testing.""" + # Create a simple quad mesh (2x2 grid) + points = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [1, 1, 0], + ], + dtype=np.uint8, + ) + + # Single quad cell + faces = np.array([4, 0, 1, 3, 2]) # quad with 4 vertices + + mesh = pv.PolyData(points, faces, force_float=False) + + # Add displacement fields for 3 timesteps + mesh.point_data["displacement_t0.000"] = np.array( + [ + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + [0, 0, 0], + ], + dtype=np.uint8, + ) + + mesh.point_data["displacement_t0.005"] = np.array( + [ + [1, 0, 0], + [1, 0, 0], + [1, 0, 0], + [1, 0, 0], + ], + dtype=np.uint8, + ) + + mesh.point_data["displacement_t0.010"] = np.array( + [ + [2, 0, 0], + [2, 0, 0], + [2, 0, 0], + [2, 0, 0], + ], + dtype=np.uint8, + ) + + # Add thickness as additional point data + mesh.point_data["thickness"] = np.array([1, 1, 1, 1], dtype=np.uint8) + + # Save to temporary file + with tempfile.NamedTemporaryFile(suffix=".vtp", delete=False) as f: + temp_path = f.name + + mesh.save(temp_path) + yield temp_path + + # Cleanup + Path(temp_path).unlink(missing_ok=True) + + +def test_load_vtp_file_basic(simple_vtp_file): + """Test basic VTP file loading.""" + pos_raw, mesh_connectivity, point_data_dict = load_vtp_file(simple_vtp_file) + + # Check positions shape: (timesteps, nodes, 3) + assert pos_raw.shape == (3, 4, 3), f"Expected shape (3, 4, 3), got {pos_raw.shape}" + + # Check mesh connectivity + assert len(mesh_connectivity) == 1, f"Expected 1 cell, got {len(mesh_connectivity)}" + assert len(mesh_connectivity[0]) == 4, ( + f"Expected quad with 4 vertices, got {len(mesh_connectivity[0])}" + ) + + # Check point data dict contains thickness + assert "thickness" in point_data_dict, "Thickness not found in point_data_dict" + assert point_data_dict["thickness"].shape == (4,), ( + f"Expected thickness shape (4,), got {point_data_dict['thickness'].shape}" + ) + + +def test_load_vtp_file_displacements(simple_vtp_file): + """Test that displacements are correctly applied.""" + pos_raw, _, _ = load_vtp_file(simple_vtp_file) + + # First timestep should be reference coords (displacement = 0) + expected_t0 = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [0, 1, 0], + [1, 1, 0], + ] + ) + np.testing.assert_array_almost_equal(pos_raw[0], expected_t0, decimal=5) + + # Second timestep should include displacement + expected_t1 = expected_t0 + np.array([[1, 0, 0]] * 4) + np.testing.assert_array_almost_equal(pos_raw[1], expected_t1, decimal=5) + + # Third timestep + expected_t2 = expected_t0 + np.array([[2, 0, 0]] * 4) + np.testing.assert_array_almost_equal(pos_raw[2], expected_t2, decimal=5) + + +def test_extract_mesh_connectivity(): + """Test mesh connectivity extraction from PolyData.""" + points = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [1, 1, 0], + [0, 1, 0], + ] + ) + + # Create a single quad + faces = np.array([4, 0, 1, 2, 3]) + poly = pv.PolyData(points, faces, force_float=False) + + connectivity = extract_mesh_connectivity_from_polydata(poly) + + assert len(connectivity) == 1, f"Expected 1 cell, got {len(connectivity)}" + assert len(connectivity[0]) == 4, f"Expected 4 vertices, got {len(connectivity[0])}" + assert connectivity[0] == [0, 1, 2, 3], ( + f"Expected [0, 1, 2, 3], got {connectivity[0]}" + ) + + +def test_build_edges_from_mesh_connectivity(): + """Test edge building from mesh connectivity.""" + # Single quad: should produce 4 edges + mesh_connectivity = [[0, 1, 2, 3]] + edges = build_edges_from_mesh_connectivity(mesh_connectivity) + + expected_edges = {(0, 1), (1, 2), (2, 3), (0, 3)} + assert edges == expected_edges, f"Expected {expected_edges}, got {edges}" + + +def test_point_data_extraction(simple_vtp_file): + """Test that non-displacement point data is extracted correctly.""" + _, _, point_data_dict = load_vtp_file(simple_vtp_file) + + # Should have thickness + assert "thickness" in point_data_dict, "Thickness not in point_data_dict" + + # Should NOT have displacement fields + assert "displacement_t0.000" not in point_data_dict, ( + "Displacement fields should not be in point_data_dict" + ) + assert "displacement_t0.005" not in point_data_dict, ( + "Displacement fields should not be in point_data_dict" + ) + + # Check thickness values + expected_thickness = np.array([1, 1, 1, 1], dtype=np.uint8) + np.testing.assert_array_almost_equal( + point_data_dict["thickness"], expected_thickness, decimal=5 + ) + + +def test_missing_displacement_fields(): + """Test that missing displacement fields raises appropriate error.""" + # Create VTP without displacement fields + points = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]]) + faces = np.array([3, 0, 1, 2]) + mesh = pv.PolyData(points, faces, force_float=False) + + with tempfile.NamedTemporaryFile(suffix=".vtp", delete=False) as f: + temp_path = f.name + + mesh.save(temp_path) + + try: + with pytest.raises(ValueError, match="No displacement fields found"): + load_vtp_file(temp_path) + finally: + Path(temp_path).unlink(missing_ok=True) + + +def test_empty_mesh_connectivity(): + """Test edge building with empty connectivity.""" + mesh_connectivity = [] + edges = build_edges_from_mesh_connectivity(mesh_connectivity) + + assert len(edges) == 0, f"Expected 0 edges for empty connectivity, got {len(edges)}" diff --git a/examples/structural_mechanics/crash/vtp_reader.py b/examples/structural_mechanics/crash/vtp_reader.py index 3cb3e88a06..cb1d99cb03 100644 --- a/examples/structural_mechanics/crash/vtp_reader.py +++ b/examples/structural_mechanics/crash/vtp_reader.py @@ -55,13 +55,14 @@ def extract_mesh_connectivity_from_polydata(poly: pv.PolyData): def load_vtp_file(vtp_path): - """Load positions over time and connectivity from a single VTP file. + """Load positions over time, connectivity, and other point data from a single VTP file. Expects displacement fields in point_data named like: - displacement_t0.000, displacement_t0.005, ..., displacement_t0.100 Returns: pos_raw: (timesteps, num_nodes, 3) absolute positions (coords + displacement_t) mesh_connectivity: list[list[int]] + point_data_dict: dict of other point data arrays (e.g., thickness) """ poly = pv.read(vtp_path) if not isinstance(poly, pv.PolyData): @@ -104,7 +105,14 @@ def natural_key(name): pos_raw = np.stack(pos_list, axis=0) mesh_connectivity = extract_mesh_connectivity_from_polydata(poly) - return pos_raw, mesh_connectivity + + # Extract all other point data fields (not displacement fields) + point_data_dict = {} + for name in poly.point_data.keys(): + if not name.startswith("displacement_"): + point_data_dict[name] = np.asarray(poly.point_data[name]) + + return pos_raw, mesh_connectivity, point_data_dict def build_edges_from_mesh_connectivity(mesh_connectivity): @@ -178,7 +186,7 @@ def process_vtp_data(data_dir, num_samples=2, write_vtp=False, logger=None): output_dir = f"./output_{os.path.splitext(os.path.basename(vtp_path))[0]}" os.makedirs(output_dir, exist_ok=True) - pos_raw, mesh_connectivity = load_vtp_file(vtp_path) + pos_raw, mesh_connectivity, point_data_dict = load_vtp_file(vtp_path) # Use unfiltered data filtered_pos_raw = pos_raw @@ -200,7 +208,11 @@ def process_vtp_data(data_dir, num_samples=2, write_vtp=False, logger=None): write_vtp=write_vtp, logger=logger, ) - point_data_all.append({"coords": mesh_pos_all}) + + # Create record with coords and all other point data fields + record = {"coords": mesh_pos_all} + record.update(point_data_dict) # Add thickness and any other fields + point_data_all.append(record) processed_runs += 1 if processed_runs >= num_samples: