diff --git a/docs/figures/density_sensitivity.png b/docs/figures/density_sensitivity.png new file mode 100644 index 0000000..992e2fa Binary files /dev/null and b/docs/figures/density_sensitivity.png differ diff --git a/docs/figures/modal_frequencies.png b/docs/figures/modal_frequencies.png new file mode 100644 index 0000000..fa0d3f8 Binary files /dev/null and b/docs/figures/modal_frequencies.png differ diff --git a/docs/figures/section_height_sensitivity.png b/docs/figures/section_height_sensitivity.png new file mode 100644 index 0000000..ce437e8 Binary files /dev/null and b/docs/figures/section_height_sensitivity.png differ diff --git a/docs/figures/thickness_sensitivity.png b/docs/figures/thickness_sensitivity.png new file mode 100644 index 0000000..0d5036a Binary files /dev/null and b/docs/figures/thickness_sensitivity.png differ diff --git a/docs/figures/youngs_modulus_sensitivity.png b/docs/figures/youngs_modulus_sensitivity.png new file mode 100644 index 0000000..02e2524 Binary files /dev/null and b/docs/figures/youngs_modulus_sensitivity.png differ diff --git a/docs/validation/density_sensitivity.csv b/docs/validation/density_sensitivity.csv new file mode 100644 index 0000000..1600635 --- /dev/null +++ b/docs/validation/density_sensitivity.csv @@ -0,0 +1,6 @@ +parameter_value,mode_1_frequency_hz +6280.0,55.23184 +7065.0,52.07308 +7850.0,49.40086 +8635.0,47.10187 +9420.0,45.09661 diff --git a/docs/validation/modal_frequencies.csv b/docs/validation/modal_frequencies.csv new file mode 100644 index 0000000..bacfe59 --- /dev/null +++ b/docs/validation/modal_frequencies.csv @@ -0,0 +1,11 @@ +mode,frequency_hz +1,49.40086 +2,84.15125 +3,305.2293 +4,506.5651 +5,692.4157 +6,836.315 +7,1292.885 +8,1339.754 +9,1591.106 +10,2078.087 diff --git a/docs/validation/section_height_sensitivity.csv b/docs/validation/section_height_sensitivity.csv new file mode 100644 index 0000000..259836a --- /dev/null +++ b/docs/validation/section_height_sensitivity.csv @@ -0,0 +1,6 @@ +parameter_value,mode_1_frequency_hz +0.032,40.01148 +0.036000000000000004,44.74873 +0.04,49.40086 +0.044000000000000004,53.97453 +0.048,58.47559 diff --git a/docs/validation/thickness_sensitivity.csv b/docs/validation/thickness_sensitivity.csv new file mode 100644 index 0000000..6c3d077 --- /dev/null +++ b/docs/validation/thickness_sensitivity.csv @@ -0,0 +1,6 @@ +parameter_value,mode_1_frequency_hz +0.0012000000000000001,49.76519 +0.00135,49.58304 +0.0015,49.40086 +0.0016500000000000002,49.21864 +0.0018,49.03639 diff --git a/docs/validation/youngs_modulus_sensitivity.csv b/docs/validation/youngs_modulus_sensitivity.csv new file mode 100644 index 0000000..6ccb478 --- /dev/null +++ b/docs/validation/youngs_modulus_sensitivity.csv @@ -0,0 +1,6 @@ +parameter_value,mode_1_frequency_hz +168000000000.0,44.18547 +189000000000.0,46.86577 +210000000000.0,49.40086 +231000000000.00003,51.81206 +252000000000.0,54.11593 diff --git a/scripts/run_modal_analysis.py b/scripts/run_modal_analysis.py new file mode 100644 index 0000000..48399a6 --- /dev/null +++ b/scripts/run_modal_analysis.py @@ -0,0 +1,113 @@ +import csv +from pathlib import Path + +import matplotlib.pyplot as plt + +from bodysimpy.analysis.modal_validation import ( + validate_modal_response, +) +from bodysimpy.config.loader import load_config +from bodysimpy.modeling.crossmember import ( + build_crossmember_model, +) + + +def main() -> None: + config = load_config("configs/baseline_crossmember.yaml") + + model = build_crossmember_model(config) + + result = validate_modal_response( + model, + modes=10, + ) + + table_path = Path("docs/validation/modal_frequencies.csv") + + table_path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + with table_path.open( + "w", + encoding="utf-8", + newline="", + ) as stream: + writer = csv.writer(stream) + + writer.writerow( + [ + "mode", + "frequency_hz", + ] + ) + + for mode, frequency in enumerate( + result.fea_frequencies_hz, + start=1, + ): + writer.writerow( + [ + mode, + frequency, + ] + ) + + modes = list( + range( + 1, + len(result.fea_frequencies_hz) + 1, + ) + ) + + figure_path = Path("docs/figures/modal_frequencies.png") + + figure_path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + plt.figure(figsize=(8, 5)) + + plt.plot( + modes, + result.fea_frequencies_hz, + marker="o", + ) + + plt.xlabel("Mode number") + plt.ylabel("Natural frequency [Hz]") + plt.title("BodySimPy Baseline Natural Frequencies") + + plt.xticks(modes) + plt.grid(True) + plt.tight_layout() + + plt.savefig( + figure_path, + dpi=200, + ) + + plt.close() + + print() + print("BodySimPy Modal Analysis") + print("-" * 55) + + for mode, frequency in enumerate( + result.fea_frequencies_hz, + start=1, + ): + print(f"Mode {mode:>2}: {frequency:>12.4f} Hz") + + print() + print(f"Analytical mode 1: {result.analytical_mode_1_hz:.4f} Hz") + + print(f"FEA mode 1: {result.fea_mode_1_hz:.4f} Hz") + + print(f"Mode-1 error: {result.mode_1_error_percent:.4f} %") + + +if __name__ == "__main__": + main() diff --git a/scripts/run_modal_sensitivity.py b/scripts/run_modal_sensitivity.py new file mode 100644 index 0000000..1146431 --- /dev/null +++ b/scripts/run_modal_sensitivity.py @@ -0,0 +1,145 @@ +import csv +from pathlib import Path + +import matplotlib.pyplot as plt + +from bodysimpy.analysis.modal_sensitivity import ( + ModalSensitivityPoint, + density_sensitivity, + section_height_sensitivity, + thickness_sensitivity, + youngs_modulus_sensitivity, +) +from bodysimpy.config.loader import load_config +from bodysimpy.modeling.crossmember import ( + build_crossmember_model, +) + +MULTIPLIERS = ( + 0.8, + 0.9, + 1.0, + 1.1, + 1.2, +) + + +def write_study( + *, + name: str, + x_label: str, + points: tuple[ModalSensitivityPoint, ...], + scale: float, +) -> None: + csv_path = Path(f"docs/validation/{name}_sensitivity.csv") + + csv_path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + with csv_path.open( + "w", + encoding="utf-8", + newline="", + ) as stream: + writer = csv.writer(stream) + + writer.writerow( + [ + "parameter_value", + "mode_1_frequency_hz", + ] + ) + + for point in points: + writer.writerow( + [ + point.parameter_value, + point.mode_1_frequency_hz, + ] + ) + + x_values = [point.parameter_value * scale for point in points] + + frequencies = [point.mode_1_frequency_hz for point in points] + + figure_path = Path(f"docs/figures/{name}_sensitivity.png") + + figure_path.parent.mkdir( + parents=True, + exist_ok=True, + ) + + plt.figure(figsize=(7, 5)) + + plt.plot( + x_values, + frequencies, + marker="o", + ) + + plt.xlabel(x_label) + plt.ylabel("Mode-1 natural frequency [Hz]") + + plt.title(f"Mode-1 Sensitivity — {name.replace('_', ' ').title()}") + + plt.grid(True) + plt.tight_layout() + + plt.savefig( + figure_path, + dpi=200, + ) + + plt.close() + + +def main() -> None: + config = load_config("configs/baseline_crossmember.yaml") + + model = build_crossmember_model(config) + + write_study( + name="thickness", + x_label="Wall thickness [mm]", + points=thickness_sensitivity( + model, + multipliers=MULTIPLIERS, + ), + scale=1000.0, + ) + + write_study( + name="youngs_modulus", + x_label="Young's modulus [GPa]", + points=youngs_modulus_sensitivity( + model, + multipliers=MULTIPLIERS, + ), + scale=1e-9, + ) + + write_study( + name="density", + x_label="Density [kg/m³]", + points=density_sensitivity( + model, + multipliers=MULTIPLIERS, + ), + scale=1.0, + ) + + write_study( + name="section_height", + x_label="Section height [mm]", + points=section_height_sensitivity( + model, + multipliers=MULTIPLIERS, + ), + scale=1000.0, + ) + + +if __name__ == "__main__": + main() diff --git a/src/bodysimpy/analysis/modal_sensitivity.py b/src/bodysimpy/analysis/modal_sensitivity.py new file mode 100644 index 0000000..c96c400 --- /dev/null +++ b/src/bodysimpy/analysis/modal_sensitivity.py @@ -0,0 +1,149 @@ +from dataclasses import dataclass, replace + +from bodysimpy.domain.structural_model import ( + StructuralModel, +) +from bodysimpy.solvers.calculix import ( + CalculiXSolver, +) + + +@dataclass(frozen=True, slots=True) +class ModalSensitivityPoint: + parameter_value: float + mode_1_frequency_hz: float + + +def _solve_mode_1( + model: StructuralModel, +) -> float: + result = CalculiXSolver().run_modal( + model, + modes=1, + ) + + return result.natural_frequencies_hz[0] + + +def thickness_sensitivity( + model: StructuralModel, + *, + multipliers: tuple[float, ...], +) -> tuple[ModalSensitivityPoint, ...]: + points: list[ModalSensitivityPoint] = [] + + for multiplier in multipliers: + thickness = model.section.thickness_m * multiplier + + section = replace( + model.section, + thickness_m=thickness, + ) + + variant = replace( + model, + name=(f"{model.name}_thickness_{multiplier:.2f}"), + section=section, + ) + + points.append( + ModalSensitivityPoint( + parameter_value=thickness, + mode_1_frequency_hz=(_solve_mode_1(variant)), + ) + ) + + return tuple(points) + + +def youngs_modulus_sensitivity( + model: StructuralModel, + *, + multipliers: tuple[float, ...], +) -> tuple[ModalSensitivityPoint, ...]: + points: list[ModalSensitivityPoint] = [] + + for multiplier in multipliers: + youngs_modulus = model.material.youngs_modulus_pa * multiplier + + material = replace( + model.material, + youngs_modulus_pa=youngs_modulus, + ) + + variant = replace( + model, + name=(f"{model.name}_youngs_{multiplier:.2f}"), + material=material, + ) + + points.append( + ModalSensitivityPoint( + parameter_value=youngs_modulus, + mode_1_frequency_hz=(_solve_mode_1(variant)), + ) + ) + + return tuple(points) + + +def density_sensitivity( + model: StructuralModel, + *, + multipliers: tuple[float, ...], +) -> tuple[ModalSensitivityPoint, ...]: + points: list[ModalSensitivityPoint] = [] + + for multiplier in multipliers: + density = model.material.density_kg_m3 * multiplier + + material = replace( + model.material, + density_kg_m3=density, + ) + + variant = replace( + model, + name=(f"{model.name}_density_{multiplier:.2f}"), + material=material, + ) + + points.append( + ModalSensitivityPoint( + parameter_value=density, + mode_1_frequency_hz=(_solve_mode_1(variant)), + ) + ) + + return tuple(points) + + +def section_height_sensitivity( + model: StructuralModel, + *, + multipliers: tuple[float, ...], +) -> tuple[ModalSensitivityPoint, ...]: + points: list[ModalSensitivityPoint] = [] + + for multiplier in multipliers: + height = model.section.height_m * multiplier + + section = replace( + model.section, + height_m=height, + ) + + variant = replace( + model, + name=(f"{model.name}_height_{multiplier:.2f}"), + section=section, + ) + + points.append( + ModalSensitivityPoint( + parameter_value=height, + mode_1_frequency_hz=(_solve_mode_1(variant)), + ) + ) + + return tuple(points) diff --git a/src/bodysimpy/analysis/modal_validation.py b/src/bodysimpy/analysis/modal_validation.py new file mode 100644 index 0000000..3323c2f --- /dev/null +++ b/src/bodysimpy/analysis/modal_validation.py @@ -0,0 +1,56 @@ +from dataclasses import dataclass + +from bodysimpy.analysis.validation import ( + relative_error_percent, +) +from bodysimpy.domain.structural_model import ( + StructuralModel, +) +from bodysimpy.solvers.analytical import ( + AnalyticalSolver, +) +from bodysimpy.solvers.calculix import ( + CalculiXSolver, +) + + +@dataclass(frozen=True, slots=True) +class ModalValidationResult: + analytical_mode_1_hz: float + fea_mode_1_hz: float + mode_1_error_percent: float + fea_frequencies_hz: tuple[float, ...] + + +def validate_modal_response( + model: StructuralModel, + *, + modes: int = 10, +) -> ModalValidationResult: + """Validate first FE natural frequency against beam theory.""" + + analytical = AnalyticalSolver().run(model) + + if not analytical.natural_frequencies_hz: + raise ValueError("Analytical solver returned no natural frequency.") + + modal = CalculiXSolver().run_modal( + model, + modes=modes, + ) + + analytical_mode_1 = analytical.natural_frequencies_hz[0] + + fea_mode_1 = modal.natural_frequencies_hz[0] + + error = relative_error_percent( + reference=analytical_mode_1, + value=fea_mode_1, + ) + + return ModalValidationResult( + analytical_mode_1_hz=analytical_mode_1, + fea_mode_1_hz=fea_mode_1, + mode_1_error_percent=error, + fea_frequencies_hz=(modal.natural_frequencies_hz), + ) diff --git a/src/bodysimpy/domain/results.py b/src/bodysimpy/domain/results.py index 4f05492..5eccb27 100644 --- a/src/bodysimpy/domain/results.py +++ b/src/bodysimpy/domain/results.py @@ -11,3 +11,12 @@ class SimulationResult: max_axial_stress_pa: float | None = None natural_frequencies_hz: tuple[float, ...] = () work_directory: Path | None = None + + +@dataclass(frozen=True, slots=True) +class ModalResult: + """Solver-independent modal-analysis result.""" + + solver_name: str + natural_frequencies_hz: tuple[float, ...] + work_directory: Path | None = None diff --git a/src/bodysimpy/solvers/base.py b/src/bodysimpy/solvers/base.py index 437140f..340056b 100644 --- a/src/bodysimpy/solvers/base.py +++ b/src/bodysimpy/solvers/base.py @@ -1,6 +1,6 @@ from typing import Protocol -from bodysimpy.domain.results import SimulationResult +from bodysimpy.domain.results import ModalResult, SimulationResult from bodysimpy.domain.structural_model import StructuralModel @@ -10,3 +10,16 @@ class StructuralSolver(Protocol): def run(self, model: StructuralModel) -> SimulationResult: """Run a structural simulation.""" ... + + +class ModalSolver(Protocol): + """Interface implemented by modal-analysis solvers.""" + + def run_modal( + self, + model: StructuralModel, + *, + modes: int, + ) -> ModalResult: + """Run a modal analysis.""" + ... diff --git a/src/bodysimpy/solvers/calculix.py b/src/bodysimpy/solvers/calculix.py index 2669463..89fb6b4 100644 --- a/src/bodysimpy/solvers/calculix.py +++ b/src/bodysimpy/solvers/calculix.py @@ -5,18 +5,25 @@ from dataclasses import dataclass from pathlib import Path -from bodysimpy.domain.results import SimulationResult +from bodysimpy.domain.results import ModalResult, SimulationResult from bodysimpy.domain.structural_model import StructuralModel -from bodysimpy.solvers.parsers.dat_parser import parse_tip_displacement -from bodysimpy.solvers.parsers.frd_parser import parse_max_abs_stress_component +from bodysimpy.solvers.parsers.dat_parser import ( + parse_eigenfrequencies, + parse_tip_displacement, +) +from bodysimpy.solvers.parsers.frd_parser import ( + parse_max_abs_stress_component, +) class CalculiXError(RuntimeError): """Raised when a CalculiX analysis cannot be completed.""" -def build_input_deck(model: StructuralModel) -> str: - """Generate a CalculiX input deck for a cantilever box beam.""" +def _build_model_definition( + model: StructuralModel, +) -> list[str]: + """Generate the solver model definition shared by all analysis steps.""" element_count = model.mesh_elements node_count = 2 * element_count + 1 @@ -29,6 +36,7 @@ def build_input_deck(model: StructuralModel) -> str: for node_index in range(node_count): node_id = node_index + 1 + x_coordinate = model.length_m * node_index / (2 * element_count) lines.append(f"{node_id},{x_coordinate:.12e},0.0,0.0") @@ -58,7 +66,7 @@ def build_input_deck(model: StructuralModel) -> str: (f"{model.material.youngs_modulus_pa:.12e},{model.material.poisson_ratio:.12e}"), "*DENSITY", f"{model.material.density_kg_m3:.12e}", - "*BEAM SECTION,ELSET=EALL,MATERIAL=STEEL,SECTION=BOX", + ("*BEAM SECTION,ELSET=EALL,MATERIAL=STEEL,SECTION=BOX"), ( f"{model.section.height_m:.12e}," f"{model.section.width_m:.12e}," @@ -68,6 +76,21 @@ def build_input_deck(model: StructuralModel) -> str: f"{model.section.thickness_m:.12e}" ), "0.0,0.0,1.0", + ] + ) + + return lines + + +def build_input_deck( + model: StructuralModel, +) -> str: + """Generate a CalculiX static input deck.""" + + lines = _build_model_definition(model) + + lines.extend( + [ "*STEP", "*STATIC", "*CLOAD", @@ -83,6 +106,32 @@ def build_input_deck(model: StructuralModel) -> str: return "\n".join(lines) + "\n" +def build_modal_input_deck( + model: StructuralModel, + *, + modes: int, +) -> str: + """Generate a CalculiX free-vibration modal input deck.""" + + if modes <= 0: + raise ValueError("Number of requested modes must be positive.") + + lines = _build_model_definition(model) + + lines.extend( + [ + "*STEP", + "*FREQUENCY", + str(modes), + "*NODE FILE", + "U", + "*END STEP", + ] + ) + + return "\n".join(lines) + "\n" + + @dataclass(frozen=True, slots=True) class CalculiXSolver: """CalculiX CrunchiX structural solver adapter.""" @@ -91,14 +140,23 @@ class CalculiXSolver: work_root: Path = Path("results/raw") timeout_seconds: float = 60.0 - def run(self, model: StructuralModel) -> SimulationResult: + def run( + self, + model: StructuralModel, + ) -> SimulationResult: + """Run a CalculiX static structural analysis.""" + executable_path = shutil.which(self.executable) if executable_path is None: raise CalculiXError(f"CalculiX executable '{self.executable}' was not found.") work_directory = self.work_root / model.name - work_directory.mkdir(parents=True, exist_ok=True) + + work_directory.mkdir( + parents=True, + exist_ok=True, + ) job_name = "static" @@ -111,9 +169,13 @@ def run(self, model: StructuralModel) -> SimulationResult: for suffix in generated_suffixes: output_path = work_directory / f"{job_name}{suffix}" - output_path.unlink(missing_ok=True) + + output_path.unlink( + missing_ok=True, + ) input_path = work_directory / f"{job_name}.inp" + input_path.write_text( build_input_deck(model), encoding="utf-8", @@ -121,13 +183,17 @@ def run(self, model: StructuralModel) -> SimulationResult: try: completed = subprocess.run( - [executable_path, job_name], + [ + executable_path, + job_name, + ], cwd=work_directory, capture_output=True, text=True, check=False, timeout=self.timeout_seconds, ) + except subprocess.TimeoutExpired as error: raise CalculiXError("CalculiX execution exceeded the configured timeout.") from error @@ -139,6 +205,7 @@ def run(self, model: StructuralModel) -> SimulationResult: ) dat_path = work_directory / f"{job_name}.dat" + frd_path = work_directory / f"{job_name}.frd" if not dat_path.exists(): @@ -166,3 +233,100 @@ def run(self, model: StructuralModel) -> SimulationResult: max_axial_stress_pa=max_axial_stress, work_directory=work_directory, ) + + def run_modal( + self, + model: StructuralModel, + *, + modes: int, + ) -> ModalResult: + """Run a CalculiX free-vibration frequency analysis.""" + + if modes <= 0: + raise ValueError("Number of requested modes must be positive.") + + executable_path = shutil.which(self.executable) + + if executable_path is None: + raise CalculiXError(f"CalculiX executable '{self.executable}' was not found.") + + work_directory = self.work_root / f"{model.name}_modal" + + work_directory.mkdir( + parents=True, + exist_ok=True, + ) + + job_name = "modal" + + generated_suffixes = ( + ".dat", + ".frd", + ".eig", + ".sta", + ".cvg", + ) + + for suffix in generated_suffixes: + output_path = work_directory / f"{job_name}{suffix}" + + output_path.unlink( + missing_ok=True, + ) + + input_path = work_directory / f"{job_name}.inp" + + input_path.write_text( + build_modal_input_deck( + model, + modes=modes, + ), + encoding="utf-8", + ) + + try: + completed = subprocess.run( + [ + executable_path, + job_name, + ], + cwd=work_directory, + capture_output=True, + text=True, + check=False, + timeout=self.timeout_seconds, + ) + + except subprocess.TimeoutExpired as error: + raise CalculiXError( + "CalculiX modal execution exceeded the configured timeout." + ) from error + + if completed.returncode != 0: + raise CalculiXError( + "CalculiX modal execution failed.\n" + f"stdout:\n{completed.stdout}\n" + f"stderr:\n{completed.stderr}" + ) + + dat_path = work_directory / f"{job_name}.dat" + + if not dat_path.exists(): + raise CalculiXError( + "CalculiX modal analysis completed without producing the expected .dat file." + ) + + frequencies = parse_eigenfrequencies(dat_path) + + if len(frequencies) < modes: + raise CalculiXError( + "CalculiX returned fewer eigenfrequencies " + f"than requested: " + f"{len(frequencies)} < {modes}." + ) + + return ModalResult( + solver_name="calculix", + natural_frequencies_hz=(frequencies[:modes]), + work_directory=work_directory, + ) diff --git a/src/bodysimpy/solvers/parsers/dat_parser.py b/src/bodysimpy/solvers/parsers/dat_parser.py index a7ce3f8..2a6ea0d 100644 --- a/src/bodysimpy/solvers/parsers/dat_parser.py +++ b/src/bodysimpy/solvers/parsers/dat_parser.py @@ -41,3 +41,61 @@ def parse_tip_displacement( raise ValueError(f"No displacement result found for node {tip_node_id}.") return displacement + + +def parse_eigenfrequencies( + path: str | Path, +) -> tuple[float, ...]: + """Extract real eigenfrequencies in cycles per unit time from a CalculiX DAT file.""" + + dat_path = Path(path) + + frequencies: list[float] = [] + in_eigenvalue_output = False + + for line in dat_path.read_text(encoding="utf-8").splitlines(): + normalized_line = "".join(line.upper().split()) + + if "EIGENVALUEOUTPUT" in normalized_line: + in_eigenvalue_output = True + continue + + if not in_eigenvalue_output: + continue + + columns = line.split() + + if not columns: + if frequencies: + break + + continue + + if len(columns) < 4: + continue + + try: + mode_number = int(columns[0]) + + frequency_hz = float( + columns[3].replace( + "D", + "E", + ) + ) + + except ValueError: + continue + + if mode_number <= 0: + continue + + if frequency_hz <= 0.0: + continue + + frequencies.append(frequency_hz) + + if not frequencies: + raise ValueError("No eigenfrequencies were found in the CalculiX DAT file.") + + return tuple(frequencies) diff --git a/tests/integration/test_calculix_modal_solver.py b/tests/integration/test_calculix_modal_solver.py new file mode 100644 index 0000000..2297996 --- /dev/null +++ b/tests/integration/test_calculix_modal_solver.py @@ -0,0 +1,38 @@ +import shutil + +import pytest + +from bodysimpy.config.loader import load_config +from bodysimpy.modeling.crossmember import ( + build_crossmember_model, +) +from bodysimpy.solvers.calculix import ( + CalculiXSolver, +) + + +@pytest.mark.skipif( + shutil.which("ccx") is None, + reason="CalculiX is not installed.", +) +def test_calculix_modal_analysis() -> None: + config = load_config("configs/baseline_crossmember.yaml") + + model = build_crossmember_model(config) + + result = CalculiXSolver().run_modal( + model, + modes=3, + ) + + assert result.solver_name == "calculix" + + assert len(result.natural_frequencies_hz) == 3 + + assert all(frequency > 0.0 for frequency in result.natural_frequencies_hz) + + assert ( + result.natural_frequencies_hz[0] + < result.natural_frequencies_hz[1] + < result.natural_frequencies_hz[2] + ) diff --git a/tests/unit/test_calculix.py b/tests/unit/test_calculix.py index 48cb269..eca11e0 100644 --- a/tests/unit/test_calculix.py +++ b/tests/unit/test_calculix.py @@ -1,7 +1,10 @@ from bodysimpy.domain.materials import IsotropicMaterial from bodysimpy.domain.sections import RectangularHollowSection from bodysimpy.domain.structural_model import StructuralModel -from bodysimpy.solvers.calculix import build_input_deck +from bodysimpy.solvers.calculix import ( + build_input_deck, + build_modal_input_deck, +) def test_build_input_deck() -> None: @@ -33,3 +36,36 @@ def test_build_input_deck() -> None: assert "*NODE PRINT,NSET=TIP" in deck assert "*EL FILE,OUTPUT=3D" in deck assert "S,NOE" in deck + + +def test_build_modal_input_deck() -> None: + model = StructuralModel( + name="baseline_modal", + section=RectangularHollowSection( + width_m=0.080, + height_m=0.040, + thickness_m=0.0015, + ), + material=IsotropicMaterial( + name="Generic structural steel", + youngs_modulus_pa=210e9, + poisson_ratio=0.30, + density_kg_m3=7850.0, + yield_strength_pa=350e6, + ), + length_m=1.0, + tip_force_n=1000.0, + mesh_elements=20, + ) + + deck = build_modal_input_deck( + model, + modes=10, + ) + + assert "*FREQUENCY" in deck + assert "\n10\n" in deck + assert "*NODE FILE" in deck + assert "\nU\n" in deck + assert "*CLOAD" not in deck + assert "*STATIC" not in deck diff --git a/tests/unit/test_modal_dat_parser.py b/tests/unit/test_modal_dat_parser.py new file mode 100644 index 0000000..ff5ee6a --- /dev/null +++ b/tests/unit/test_modal_dat_parser.py @@ -0,0 +1,36 @@ +from pathlib import Path + +import pytest + +from bodysimpy.solvers.parsers.dat_parser import ( + parse_eigenfrequencies, +) + + +def test_parse_eigenfrequencies(tmp_path: Path) -> None: + dat_path = tmp_path / "modal.dat" + + dat_path.write_text( + """ + E I G E N V A L U E O U T P U T + + MODE NO EIGENVALUE FREQUENCY + RAD/TIME CYCLES/TIME + + 1 9.576000E+04 3.094511E+02 4.925290E+01 + 2 4.500000E+05 6.708204E+02 1.067632E+02 + 3 7.000000E+05 8.366600E+02 1.331575E+02 + +""", + encoding="utf-8", + ) + + frequencies = parse_eigenfrequencies(dat_path) + + assert frequencies == pytest.approx( + ( + 49.25290, + 106.7632, + 133.1575, + ) + )