From 1ee8b0c6a78af542e76671ebb4d16e8fd98d0721 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 7 Apr 2026 09:30:42 +0200 Subject: [PATCH 01/16] add is_loaded to beam_data --- src/py_sofistik_utils/cdb_reader/_internals/beam_data.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam_data.py b/src/py_sofistik_utils/cdb_reader/_internals/beam_data.py index 93b2b08..02224c9 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam_data.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam_data.py @@ -85,11 +85,13 @@ def __init__(self, dll: SofDll) -> None: ] ) self._dll = dll + self._is_loaded = False def clear(self) -> None: """Clear all the loaded data. """ self._data = self._data[0:0] + self._is_loaded = False def data(self, deep: bool = True) -> DataFrame: """Return the :class:`pandas.DataFrame` containing the loaded key @@ -161,6 +163,11 @@ def get( f"and quantity {quantity}!" ) from e + def is_loaded(self) -> bool: + """Return ``True`` if beam data have been loaded, ``False`` otherwise. + """ + return self._is_loaded + def load(self) -> None: """Retrieve all beam data. If the key does not exist or it is empty, a warning is raised only if ``echo_level > 0``. @@ -276,3 +283,4 @@ def load(self) -> None: # set indices for fast lookup and override existing data, if any self._data = df.set_index(["ELEM_ID"], drop=False) + self._is_loaded = True From 30c1366d83d33f7743267ac463a94e5f76aedc22 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 7 Apr 2026 09:37:29 +0200 Subject: [PATCH 02/16] add test for is_loaded in test_beam_data --- tests/cdb_reader/test_beam_data.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/cdb_reader/test_beam_data.py b/tests/cdb_reader/test_beam_data.py index 845f34f..f4ecbac 100644 --- a/tests/cdb_reader/test_beam_data.py +++ b/tests/cdb_reader/test_beam_data.py @@ -162,3 +162,11 @@ def test_get_after_clear(self) -> None: self.cdb.beams.data.load() with self.subTest(msg="Check indexes management"): self.test_get() + + def test_is_loaded(self) -> None: + with self.subTest(msg="After load"): + self. assertTrue(self.cdb.beams.data.is_loaded()) + + self.cdb.beams.data.clear() + with self.subTest(msg="After clear"): + self. assertFalse(self.cdb.beams.data.is_loaded()) From a9a78577e541a8046d83563dc3e222bea86ce865 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Thu, 30 Apr 2026 12:13:37 +0200 Subject: [PATCH 03/16] add strain energy to beam.py (WIP) --- .../cdb_reader/_internals/beam.py | 165 +++++++++++++++++- 1 file changed, 163 insertions(+), 2 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 430d0a7..af81505 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -8,6 +8,7 @@ from . beam_load import _BeamLoad from . beam_results import BeamResults from . beam_stresses import _BeamStress +from . property import _PropertyData from . sofistik_dll import SofDll @@ -31,5 +32,165 @@ def __init__(self, dll: SofDll) -> None: self.results = BeamResults(dll) self.stresses = _BeamStress(dll) - self._calculated_lc: set[int] = set() - self._data = DataFrame(columns=["LOAD_CASE", "ID", "X", "Y", "Z"]) + self._calculated_u_lc: set[int] = set() + self._dll = dll + self._u = DataFrame( + columns=[ + "ELEM_ID", + "GROUP", + "LOAD_CASE", + "U" + ] + ) + + def calculate_strain_energy(self, load_cases: int | list[int]) -> None: + """Calculate the strain energy for the given ``load_case`` numbers. + + Notes + ----- + Please refer to the documentation for the hypotheses on which the + calculation is based on. + """ + if isinstance(load_cases, int): + load_cases = [load_cases] + else: + load_cases = list(set(load_cases)) + + if not self.data.is_loaded(): + self.data.load() + + # Load properties for all unique PROP_END_1/2 + prop_nmb = ( + set(self.data.data()["PROP_END_1"].unique()) | + set(self.data.data()["PROP_END_2"].unique()) + ) + properties = _PropertyData(self._dll) + for prop in prop_nmb: + properties.load(int(prop)) + + # Load results for all load cases + self.results.load(load_cases) + + # Initialize base_data with zeros for required columns + base_data = self.results.data().reset_index(drop=True) + required_columns = ["LENGTH", "PROP_END_1", "PROP_END_2", "E", "A", "IYY", "IZZ", "U"] + base_data[required_columns] = 0.0 + + # Assign PROP_END_1/2 and LENGTH using vectorized operations + elem_ids = base_data["ELEM_ID"].unique() + base_data["PROP_END_1"] = base_data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_1") for _ in elem_ids}).fillna(1).astype(int) + base_data["PROP_END_2"] = base_data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_2") for _ in elem_ids}).fillna(1).astype(int) + base_data["LENGTH"] = base_data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) + + for prop in prop_nmb: + base_data.loc[base_data["PROP_END_1"] == prop, "A"] = properties.get_area(prop) + base_data.loc[base_data["PROP_END_1"] == prop, "IYY"] = properties.get_second_moment_of_area_yy(prop) + base_data.loc[base_data["PROP_END_1"] == prop, "IZZ"] = properties.get_second_moment_of_area_zz(prop) + base_data.loc[base_data["PROP_END_2"] == prop, "A"] = properties.get_area(prop) + base_data.loc[base_data["PROP_END_2"] == prop, "IYY"] = properties.get_second_moment_of_area_yy(prop) + base_data.loc[base_data["PROP_END_2"] == prop, "IZZ"] = properties.get_second_moment_of_area_zz(prop) + + # Calculate strain energy (vectorized) + E = 210000000 # kN/m² + base_data["U"] = (1 / (2 * E)) * ( + base_data["N"] ** 2 / base_data["A"] + + base_data["MY"] ** 2 / base_data["IYY"] + + base_data["MZ"] ** 2 / base_data["IZZ"] + ) + + # Calculate segment lengths (next_POS - current_POS) + base_data = base_data.sort_values(["LOAD_CASE", "ELEM_ID", "POS"]) + base_data["SEGMENT_LENGTH"] = base_data.groupby(["LOAD_CASE", "ELEM_ID"])["POS"].diff(-1).abs().fillna(0) + + # Calculate U for the next section (shifted U values) + base_data["U_NEXT"] = base_data.groupby(["LOAD_CASE", "ELEM_ID"])["U"].shift(-1).fillna(0) + + # Calculate strain energy per segment: (U1 + U2) * L_seg / 2 + base_data["U_SEGMENT"] = (base_data["U"] + base_data["U_NEXT"]) * base_data["SEGMENT_LENGTH"] / 2 + + # Sum strain energy across all segments for each (LOAD_CASE, ELEM_ID) + result = base_data.groupby(["LOAD_CASE", "ELEM_ID"]).agg({ + "U_SEGMENT": "sum", + "GROUP": "first" # Assuming GROUP is the same for all segments of an ELEM_ID + }).rename(columns={"U_SEGMENT": "U"}).reset_index() + + # Reorder columns to ['ELEM_ID', 'GROUP', 'LOAD_CASE', 'U'] + result = result[["ELEM_ID", "GROUP", "LOAD_CASE", "U"]] + + # Set multi-index for fast lookups + self._u = result.set_index( + ["LOAD_CASE", "ELEM_ID"], + drop=False + ).sort_index() + + self._calculated_u_lc.update(load_cases) + + def clear_strain_energy(self, load_case: int) -> None: + """Clear the calculated strain energy for the given ``load_case`` + number. + """ + if load_case not in self._calculated_u_lc: + return + + self._data = self._u[ + self._u.index.get_level_values("LOAD_CASE") != load_case + ] + self._calculated_u_lc.remove(load_case) + + def clear_all_strain_energy(self) -> None: + """Clear the calculated strain energy for all the load cases. + """ + self._u = self._u[0:0] + self._calculated_u_lc.clear() + + def get_element_strain_energy( + self, + load_case: int, + element_id: int, + default: float | None = None + ) -> float: + """Return the calculated strain energy for the given ``element_id`` and + ``load_case``. + """ + try: + return self._u.at[(load_case, element_id), "U"] # type: ignore + except (KeyError, ValueError) as e: + if default is not None: + return default + raise LookupError( + f"Strain energy entry not found for element id {element_id} " + f"and load case {load_case}!" + ) from e + + def get_group_strain_energy( + self, + load_case: int, + group: int, + default: float | None = None + ) -> float: + """Return the sum of the calculated strain energy considering all the + elements in the given ``group`` number and ``load_case``. + """ + try: + return self._u[self._u.GROUP == group].loc[load_case, "U"].sum() # type: ignore + except (KeyError, ValueError) as e: + if default is not None: + return default + raise LookupError( + f"Strain energy could not be returned for group {group} " + f"and load case {load_case}!" + ) from e + + def get_strain_energy(self, deep: bool = True) -> DataFrame: + """Return the calculated strain energy for all the loaded + ``load_case``. + + Parameters + ---------- + deep : bool, default True + When ``deep=True``, a new object will be created with a copy of the + calling object's data and indices. Modifications to the data or + indices of the copy will not be reflected in the original object + (refer to :meth:`pandas.DataFrame.copy` documentation for details). + """ + return self._u.copy(deep=deep) From 8591b831849c9e707d48fb2e18350e8d472a3c3f Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Thu, 30 Apr 2026 12:14:06 +0200 Subject: [PATCH 04/16] add test suites for strain energy --- tests/cdb_reader/test_beam_strain_energy.py | 201 ++++++++++++++++++++ 1 file changed, 201 insertions(+) create mode 100644 tests/cdb_reader/test_beam_strain_energy.py diff --git a/tests/cdb_reader/test_beam_strain_energy.py b/tests/cdb_reader/test_beam_strain_energy.py new file mode 100644 index 0000000..da9f314 --- /dev/null +++ b/tests/cdb_reader/test_beam_strain_energy.py @@ -0,0 +1,201 @@ +# standard library imports +from os import environ +from unittest import skipUnless, TestCase + +# third party library imports +from pandas import DataFrame +from pandas.testing import assert_frame_equal + +# local library specific imports +from py_sofistik_utils import SOFiSTiKCDBReader + + +CDB_PATH = environ.get("SOFISTIK_CDB_PATH") +DLL_PATH = environ.get("SOFISTIK_DLL_PATH") +VERSION = environ.get("SOFISTIK_VERSION") + +_COLUMNS = ["ELEM_ID", "GROUP", "LOAD_CASE", "U"] +_DATA = [ + [101, 10, 1000, 0.34917488541448244], + [102, 10, 1000, 1.2396276374589636], + [101, 10, 1001, 0.02940052665116039], + [102, 10, 1001, 0.0], + [101, 10, 1002, 0.3785754120656429], + [102, 10, 1002, 1.2396276374589636] +] +_DATA_DIV = [ + [101, 10, 1000, 0.3846355348887977], + [102, 10, 1000, 0.8608525260131692], + [101, 10, 1001, 0.02940052665116039], + [102, 10, 1001, 0.0], + [101, 10, 1002, 0.41403606153995803], + [102, 10, 1002, 0.8608525260131692] +] + + +@skipUnless( + all([CDB_PATH, DLL_PATH, VERSION]), + "SOFiSTiK environment variables not set!" +) +class SOFiSTiKCDBReaderBeamStrainEnergyTestSuite(TestCase): + def setUp(self) -> None: + self.cdb = SOFiSTiKCDBReader( + CDB_PATH, # type: ignore + "BEAM_STRAIN_ENERGY", + DLL_PATH, # type: ignore + int(VERSION) # type: ignore + ) + self.cdb.initialize() + self.cdb.set_echo_level(1) + self.cdb.beams.calculate_strain_energy([1000, 1001, 1002]) + + def tearDown(self) -> None: + self.cdb.close() + + def test_get_strain_energy(self) -> None: + data = DataFrame( + _DATA, + columns=_COLUMNS + ).set_index(["LOAD_CASE", "ELEM_ID"], drop=False).sort_index() + + # NOTE: + # Float values loaded from the CDB contain inherent numerical noise. + # The chosen tolerance is stricter than pandas default and reflects the + # maximum relative error observed in practice, ensuring stable and + # reproducible comparisons. + assert_frame_equal( + data, + self.cdb.beams.get_strain_energy(), + rtol=1E-7 + ) + + def test_get_element_strain_energy(self) -> None: + with self.subTest(msg="Existing entry"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1002, 101), + 0.3785754120656429 + ) + + with self.subTest(msg="Non existing entry with default"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1002, 103, -3), + -3 + ) + + with self.subTest(msg="Non existing entry without default"): + with self.assertRaises(LookupError): + self.cdb.beams.get_element_strain_energy(1002, 103) + + def test_get_group_strain_energy(self) -> None: + with self.subTest(msg="Existing entry"): + self.assertEqual( + self.cdb.beams.get_group_strain_energy(1002, 10), + 0.3785754120656429 + 1.2396276374589636 + ) + + with self.subTest(msg="Non existing entry with default"): + self.assertEqual( + self.cdb.beams.get_group_strain_energy(1002, 11, -3), + -3 + ) + + with self.subTest(msg="Non existing entry without default"): + with self.assertRaises(LookupError): + self.cdb.beams.get_group_strain_energy(1002, 11) + + def test_get_strain_energy_after_clear_all(self) -> None: + self.cdb.beams.clear_all_strain_energy() + with self.subTest(msg="Check clear all method"): + with self.assertRaises(LookupError): + self.cdb.beams.get_element_strain_energy(1001, 101) + + self.cdb.beams.calculate_strain_energy(1001) + with self.subTest(msg="Check indexes management"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1001, 101), + 0.02940052665116039 + ) + + +@skipUnless( + all([CDB_PATH, DLL_PATH, VERSION]), + "SOFiSTiK environment variables not set!" +) +class SOFiSTiKCDBReaderBeamStrainEnergyWithDIVTestSuite(TestCase): + def setUp(self) -> None: + self.cdb = SOFiSTiKCDBReader( + CDB_PATH, # type: ignore + "BEAM_STRAIN_ENERGY_WITH_DIV", + DLL_PATH, # type: ignore + int(VERSION) # type: ignore + ) + self.cdb.initialize() + self.cdb.set_echo_level(1) + self.cdb.beams.calculate_strain_energy([1000, 1001, 1002]) + + def tearDown(self) -> None: + self.cdb.close() + + def test_get_strain_energy(self) -> None: + data = DataFrame( + _DATA_DIV, + columns=_COLUMNS + ).set_index(["LOAD_CASE", "ELEM_ID"], drop=False).sort_index() + + # NOTE: + # Float values loaded from the CDB contain inherent numerical noise. + # The chosen tolerance is stricter than pandas default and reflects the + # maximum relative error observed in practice, ensuring stable and + # reproducible comparisons. + assert_frame_equal( + data, + self.cdb.beams.get_strain_energy(), + rtol=1E-7 + ) + + def test_get_element_strain_energy(self) -> None: + with self.subTest(msg="Existing entry"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1002, 101), + 0.41403606153995803 + ) + + with self.subTest(msg="Non existing entry with default"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1002, 103, -3), + -3 + ) + + with self.subTest(msg="Non existing entry without default"): + with self.assertRaises(LookupError): + self.cdb.beams.get_element_strain_energy(1002, 103) + + def test_get_group_strain_energy(self) -> None: + with self.subTest(msg="Existing entry"): + self.assertEqual( + self.cdb.beams.get_group_strain_energy(1002, 10), + 0.41403606153995803 + 0.8608525260131692 + ) + + with self.subTest(msg="Non existing entry with default"): + self.assertEqual( + self.cdb.beams.get_group_strain_energy(1002, 11, -3), + -3 + ) + + with self.subTest(msg="Non existing entry without default"): + with self.assertRaises(LookupError): + self.cdb.beams.get_group_strain_energy(1002, 11) + + def test_get_strain_energy_after_clear_all(self) -> None: + self.cdb.beams.clear_all_strain_energy() + with self.subTest(msg="Check clear all method"): + with self.assertRaises(LookupError): + self.cdb.beams.get_element_strain_energy(1001, 101) + + self.cdb.beams.calculate_strain_energy(1001) + with self.subTest(msg="Check indexes management"): + self.assertEqual( + self.cdb.beams.get_element_strain_energy(1001, 101), + 0.02940052665116039 + ) From 4a9a7afe619024cd50c8a80959014a073f8e2528 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Thu, 30 Apr 2026 13:33:38 +0200 Subject: [PATCH 05/16] add tests to docs --- docs/cdb/test_setup.rst | 1 + docs/cdb/tests/beam_strain_energy.rst | 130 ++++++++++++++++++++++++++ 2 files changed, 131 insertions(+) create mode 100644 docs/cdb/tests/beam_strain_energy.rst diff --git a/docs/cdb/test_setup.rst b/docs/cdb/test_setup.rst index 5b330b5..cff7c92 100644 --- a/docs/cdb/test_setup.rst +++ b/docs/cdb/test_setup.rst @@ -89,6 +89,7 @@ the temporary environment variable approach, open an MSYS2 MINGW64 shell, naviga tests/beam_data tests/beam_result + tests/beam_strain_energy tests/cable_data tests/cable_load tests/cable_result diff --git a/docs/cdb/tests/beam_strain_energy.rst b/docs/cdb/tests/beam_strain_energy.rst new file mode 100644 index 0000000..98df696 --- /dev/null +++ b/docs/cdb/tests/beam_strain_energy.rst @@ -0,0 +1,130 @@ +BeamStrainEnergy +---------------- + +Related test suite: ``test_beam_strain_energy.py`` + +There are two test suites for beam strain energy and therefore two expected CDB +files. The first one tests usual beam element without internal subdivisions +(DIV), while the second tests beams with subdivisions. + +Part I +"""""" + +Expected CDB file name: ``BEAM_STRAIN_ENERGY.cdb`` + +Runs with: SOFiSTiK 2025 + +Version: 1 + +.. code-block:: text + + +PROG AQUA + HEAD MATERIAL AND SECTIONS + NORM EN 199X-200X + STEE NO 1 TYPE S 355 ES 210000.0 GAM 78.5 TITL 'S355' + PROF 1 TYPE CHS 300.0 10.0 MNO 1 + PROF 2 TYPE CHS 200.0 10.0 MNO 1 + END + + +PROG SOFIMSHA + HEAD GEOMETRY REV-1-SOF-2025 + SYST 3D GDIR NEGY GDIV 10 + NODE NO 1 X 00.0 Y 0.0 Z +0.0 FIX MX,PX,PY,PZ + NODE NO 2 X 05.0 Y 0.0 Z +0.0 + NODE NO 3 X 10.0 Y 0.0 Z +0.0 FIX PP + + GRP 10 TITL 'BEAM_DIV' + BEAM NO 1 NA 1 NE 2 NCS 1 DIV 1 AHIN MT + BEAM NO 2 NA 2 NE 3 NCS 2.1 EHIN NMYMZ + END + + +PROG SOFILOAD + HEAD LOADS + LC 20 TITL 'LOAD-1' + NODE 2 TYPE PYY -30.0 + LC 30 TITL 'LOAD-2' + NODE 2 TYPE PXX +150.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1000 DLY 0.0 TITL 'LC 1000' + LCC 20 FACT 1.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1001 DLY 0.0 TITL 'LC 1001' + LCC 30 FACT 1.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1002 DLY 0.0 TITL 'LC 1002' + LCC 20 FACT 1.0 + LCC 30 FACT 1.0 + END + +Part II +"""""" + +Expected CDB file name: ``BEAM_STRAIN_ENERGY_WITH_DIV.cdb`` + +Runs with: SOFiSTiK 2025 + +Version: 1 + +.. code-block:: text + + +PROG AQUA + HEAD MATERIAL AND SECTIONS + NORM EN 199X-200X + STEE NO 1 TYPE S 355 ES 210000.0 GAM 78.5 TITL 'S355' + PROF 1 TYPE CHS 300.0 10.0 MNO 1 + PROF 2 TYPE CHS 200.0 10.0 MNO 1 + END + + +PROG SOFIMSHA + HEAD GEOMETRY REV-1-SOF-2025 + SYST 3D GDIR NEGY GDIV 10 + NODE NO 1 X 00.0 Y 0.0 Z +0.0 FIX MX,PX,PY,PZ + NODE NO 2 X 05.0 Y 0.0 Z +0.0 + NODE NO 3 X 10.0 Y 0.0 Z +0.0 FIX PP + + GRP 10 TITL 'BEAM_DIV' + BEAM NO 1 NA 1 NE 2 NCS 1 DIV 5 AHIN MT + BEAM NO 2 NA 2 NE 3 NCS 2.1 EHIN NMYMZ + END + + +PROG SOFILOAD + HEAD LOADS + LC 20 TITL 'LOAD-1' + BEAM 101 TYPE PYY -10.0 + LC 30 TITL 'LOAD-2' + NODE 2 TYPE PXX +150.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1000 DLY 0.0 TITL 'LC 1000' + LCC 20 FACT 1.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1001 DLY 0.0 TITL 'LC 1001' + LCC 30 FACT 1.0 + END + + +PROG ASE + HEAD LINEAR ANALYSIS + SYST PROB LINE + LC 1002 DLY 0.0 TITL 'LC 1002' + LCC 20 FACT 1.0 + LCC 30 FACT 1.0 + END From 2c0211d622122550f355ba30bd73cc277f6af534 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Mon, 4 May 2026 11:44:29 +0200 Subject: [PATCH 06/16] typo in beam_strain_energy.rst --- docs/cdb/tests/beam_strain_energy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/cdb/tests/beam_strain_energy.rst b/docs/cdb/tests/beam_strain_energy.rst index 98df696..a338e63 100644 --- a/docs/cdb/tests/beam_strain_energy.rst +++ b/docs/cdb/tests/beam_strain_energy.rst @@ -69,7 +69,7 @@ Version: 1 END Part II -"""""" +""""""" Expected CDB file name: ``BEAM_STRAIN_ENERGY_WITH_DIV.cdb`` From 0ad383a0e9e76227d583c277bf36e6aea36e24ea Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Mon, 4 May 2026 11:49:02 +0200 Subject: [PATCH 07/16] update beam.py after #62 --- .../cdb_reader/_internals/beam.py | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index af81505..d44999f 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -8,7 +8,7 @@ from . beam_load import _BeamLoad from . beam_results import BeamResults from . beam_stresses import _BeamStress -from . property import _PropertyData +from . cross_section_data import CrossSectionalData from . sofistik_dll import SofDll @@ -60,13 +60,12 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: self.data.load() # Load properties for all unique PROP_END_1/2 - prop_nmb = ( + prop_nmb: list[int] = list( set(self.data.data()["PROP_END_1"].unique()) | set(self.data.data()["PROP_END_2"].unique()) ) - properties = _PropertyData(self._dll) - for prop in prop_nmb: - properties.load(int(prop)) + properties = CrossSectionalData(self._dll) + properties.load(prop_nmb) # Load results for all load cases self.results.load(load_cases) @@ -83,16 +82,17 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: base_data["LENGTH"] = base_data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) for prop in prop_nmb: - base_data.loc[base_data["PROP_END_1"] == prop, "A"] = properties.get_area(prop) - base_data.loc[base_data["PROP_END_1"] == prop, "IYY"] = properties.get_second_moment_of_area_yy(prop) - base_data.loc[base_data["PROP_END_1"] == prop, "IZZ"] = properties.get_second_moment_of_area_zz(prop) - base_data.loc[base_data["PROP_END_2"] == prop, "A"] = properties.get_area(prop) - base_data.loc[base_data["PROP_END_2"] == prop, "IYY"] = properties.get_second_moment_of_area_yy(prop) - base_data.loc[base_data["PROP_END_2"] == prop, "IZZ"] = properties.get_second_moment_of_area_zz(prop) + base_data.loc[base_data["PROP_END_1"] == prop, "A"] = properties.get(prop, "A") + base_data.loc[base_data["PROP_END_1"] == prop, "IYY"] = properties.get(prop, "IY") + base_data.loc[base_data["PROP_END_1"] == prop, "IZZ"] = properties.get(prop, "IZ") + base_data.loc[base_data["PROP_END_1"] == prop, "E"] = properties.get(prop, "EM") + base_data.loc[base_data["PROP_END_2"] == prop, "A"] = properties.get(prop, "A") + base_data.loc[base_data["PROP_END_2"] == prop, "IYY"] = properties.get(prop, "IY") + base_data.loc[base_data["PROP_END_2"] == prop, "IZZ"] = properties.get(prop, "IZ") + base_data.loc[base_data["PROP_END_2"] == prop, "E"] = properties.get(prop, "EM") # Calculate strain energy (vectorized) - E = 210000000 # kN/m² - base_data["U"] = (1 / (2 * E)) * ( + base_data["U"] = (1 / (2 * base_data["E"])) * ( base_data["N"] ** 2 / base_data["A"] + base_data["MY"] ** 2 / base_data["IYY"] + base_data["MZ"] ** 2 / base_data["IZZ"] From f42d6c024a09973b335470d1dc5895be5eec755e Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Mon, 4 May 2026 12:44:56 +0200 Subject: [PATCH 08/16] CI trigger From 2657697a7c5bcced05b6df8a2cf73080f8fdd5b2 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 10:19:02 +0200 Subject: [PATCH 09/16] fix wrong type in self.key_exist in beam.py --- src/py_sofistik_utils/cdb_reader/_internals/beam.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index d44999f..9e2f61d 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -65,7 +65,7 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: set(self.data.data()["PROP_END_2"].unique()) ) properties = CrossSectionalData(self._dll) - properties.load(prop_nmb) + properties.load([int(_) for _ in prop_nmb]) # Load results for all load cases self.results.load(load_cases) From 7e4ac6285705a6bbd38f5aea047622d370bc1c8a Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 11:48:47 +0200 Subject: [PATCH 10/16] make column names consistent with cross_section --- .../cdb_reader/_internals/beam.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 9e2f61d..e928736 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -72,7 +72,7 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: # Initialize base_data with zeros for required columns base_data = self.results.data().reset_index(drop=True) - required_columns = ["LENGTH", "PROP_END_1", "PROP_END_2", "E", "A", "IYY", "IZZ", "U"] + required_columns = ["LENGTH", "PROP_END_1", "PROP_END_2", "EM", "A", "IY", "IZ", "U"] base_data[required_columns] = 0.0 # Assign PROP_END_1/2 and LENGTH using vectorized operations @@ -83,19 +83,19 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: for prop in prop_nmb: base_data.loc[base_data["PROP_END_1"] == prop, "A"] = properties.get(prop, "A") - base_data.loc[base_data["PROP_END_1"] == prop, "IYY"] = properties.get(prop, "IY") - base_data.loc[base_data["PROP_END_1"] == prop, "IZZ"] = properties.get(prop, "IZ") - base_data.loc[base_data["PROP_END_1"] == prop, "E"] = properties.get(prop, "EM") + base_data.loc[base_data["PROP_END_1"] == prop, "IY"] = properties.get(prop, "IY") + base_data.loc[base_data["PROP_END_1"] == prop, "IZ"] = properties.get(prop, "IZ") + base_data.loc[base_data["PROP_END_1"] == prop, "EM"] = properties.get(prop, "EM") base_data.loc[base_data["PROP_END_2"] == prop, "A"] = properties.get(prop, "A") - base_data.loc[base_data["PROP_END_2"] == prop, "IYY"] = properties.get(prop, "IY") - base_data.loc[base_data["PROP_END_2"] == prop, "IZZ"] = properties.get(prop, "IZ") - base_data.loc[base_data["PROP_END_2"] == prop, "E"] = properties.get(prop, "EM") + base_data.loc[base_data["PROP_END_2"] == prop, "IY"] = properties.get(prop, "IY") + base_data.loc[base_data["PROP_END_2"] == prop, "IZ"] = properties.get(prop, "IZ") + base_data.loc[base_data["PROP_END_2"] == prop, "EM"] = properties.get(prop, "EM") # Calculate strain energy (vectorized) - base_data["U"] = (1 / (2 * base_data["E"])) * ( + base_data["U"] = (1 / (2 * base_data["EM"])) * ( base_data["N"] ** 2 / base_data["A"] + - base_data["MY"] ** 2 / base_data["IYY"] + - base_data["MZ"] ** 2 / base_data["IZZ"] + base_data["MY"] ** 2 / base_data["IY"] + + base_data["MZ"] ** 2 / base_data["IZ"] ) # Calculate segment lengths (next_POS - current_POS) From af9738a27ea718d9bd5c506279d8a9769e2baca1 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 11:54:26 +0200 Subject: [PATCH 11/16] clean-up --- .../cdb_reader/_internals/beam.py | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index e928736..94c3c35 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -70,48 +70,48 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: # Load results for all load cases self.results.load(load_cases) - # Initialize base_data with zeros for required columns - base_data = self.results.data().reset_index(drop=True) + # Initialize data with zeros for required columns + data = self.results.data().reset_index(drop=True) required_columns = ["LENGTH", "PROP_END_1", "PROP_END_2", "EM", "A", "IY", "IZ", "U"] - base_data[required_columns] = 0.0 + data[required_columns] = 0.0 # Assign PROP_END_1/2 and LENGTH using vectorized operations - elem_ids = base_data["ELEM_ID"].unique() - base_data["PROP_END_1"] = base_data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_1") for _ in elem_ids}).fillna(1).astype(int) - base_data["PROP_END_2"] = base_data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_2") for _ in elem_ids}).fillna(1).astype(int) - base_data["LENGTH"] = base_data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) - - for prop in prop_nmb: - base_data.loc[base_data["PROP_END_1"] == prop, "A"] = properties.get(prop, "A") - base_data.loc[base_data["PROP_END_1"] == prop, "IY"] = properties.get(prop, "IY") - base_data.loc[base_data["PROP_END_1"] == prop, "IZ"] = properties.get(prop, "IZ") - base_data.loc[base_data["PROP_END_1"] == prop, "EM"] = properties.get(prop, "EM") - base_data.loc[base_data["PROP_END_2"] == prop, "A"] = properties.get(prop, "A") - base_data.loc[base_data["PROP_END_2"] == prop, "IY"] = properties.get(prop, "IY") - base_data.loc[base_data["PROP_END_2"] == prop, "IZ"] = properties.get(prop, "IZ") - base_data.loc[base_data["PROP_END_2"] == prop, "EM"] = properties.get(prop, "EM") + elem_ids = data["ELEM_ID"].unique() + data["PROP_END_1"] = data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_1") for _ in elem_ids}).fillna(1).astype(int) + data["PROP_END_2"] = data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_2") for _ in elem_ids}).fillna(1).astype(int) + data["LENGTH"] = data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) + + for p in prop_nmb: + data.loc[data["PROP_END_1"] == p, "A"] = properties.get(p, "A") + data.loc[data["PROP_END_1"] == p, "IY"] = properties.get(p, "IY") + data.loc[data["PROP_END_1"] == p, "IZ"] = properties.get(p, "IZ") + data.loc[data["PROP_END_1"] == p, "EM"] = properties.get(p, "EM") + data.loc[data["PROP_END_2"] == p, "A"] = properties.get(p, "A") + data.loc[data["PROP_END_2"] == p, "IY"] = properties.get(p, "IY") + data.loc[data["PROP_END_2"] == p, "IZ"] = properties.get(p, "IZ") + data.loc[data["PROP_END_2"] == p, "EM"] = properties.get(p, "EM") # Calculate strain energy (vectorized) - base_data["U"] = (1 / (2 * base_data["EM"])) * ( - base_data["N"] ** 2 / base_data["A"] + - base_data["MY"] ** 2 / base_data["IY"] + - base_data["MZ"] ** 2 / base_data["IZ"] + data["U"] = (1 / (2 * data["EM"])) * ( + data["N"] ** 2 / data["A"] + + data["MY"] ** 2 / data["IY"] + + data["MZ"] ** 2 / data["IZ"] ) # Calculate segment lengths (next_POS - current_POS) - base_data = base_data.sort_values(["LOAD_CASE", "ELEM_ID", "POS"]) - base_data["SEGMENT_LENGTH"] = base_data.groupby(["LOAD_CASE", "ELEM_ID"])["POS"].diff(-1).abs().fillna(0) + data = data.sort_values(["LOAD_CASE", "ELEM_ID", "POS"]) + data["SEGMENT_LENGTH"] = data.groupby(["LOAD_CASE", "ELEM_ID"])["POS"].diff(-1).abs().fillna(0) # Calculate U for the next section (shifted U values) - base_data["U_NEXT"] = base_data.groupby(["LOAD_CASE", "ELEM_ID"])["U"].shift(-1).fillna(0) + data["U_NEXT"] = data.groupby(["LOAD_CASE", "ELEM_ID"])["U"].shift(-1).fillna(0) # Calculate strain energy per segment: (U1 + U2) * L_seg / 2 - base_data["U_SEGMENT"] = (base_data["U"] + base_data["U_NEXT"]) * base_data["SEGMENT_LENGTH"] / 2 + data["U_SEGMENT"] = (data["U"] + data["U_NEXT"]) * data["SEGMENT_LENGTH"] / 2 # Sum strain energy across all segments for each (LOAD_CASE, ELEM_ID) - result = base_data.groupby(["LOAD_CASE", "ELEM_ID"]).agg({ + result = data.groupby(["LOAD_CASE", "ELEM_ID"]).agg({ "U_SEGMENT": "sum", - "GROUP": "first" # Assuming GROUP is the same for all segments of an ELEM_ID + "GROUP": "first" }).rename(columns={"U_SEGMENT": "U"}).reset_index() # Reorder columns to ['ELEM_ID', 'GROUP', 'LOAD_CASE', 'U'] From eb9a351bdb30127bc58386ac54af74eae680b5d8 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 12:07:00 +0200 Subject: [PATCH 12/16] use pandas.eval to improve performance --- src/py_sofistik_utils/cdb_reader/_internals/beam.py | 8 +++----- tests/cdb_reader/test_beam_strain_energy.py | 12 ++++++------ 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 94c3c35..292a7d6 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -91,11 +91,9 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: data.loc[data["PROP_END_2"] == p, "IZ"] = properties.get(p, "IZ") data.loc[data["PROP_END_2"] == p, "EM"] = properties.get(p, "EM") - # Calculate strain energy (vectorized) - data["U"] = (1 / (2 * data["EM"])) * ( - data["N"] ** 2 / data["A"] + - data["MY"] ** 2 / data["IY"] + - data["MZ"] ** 2 / data["IZ"] + # Calculate strain energy + data["U"] = data.eval( + "(N ** 2 / A + MY ** 2 / IY + MZ ** 2 / IZ) / (2 * EM)" ) # Calculate segment lengths (next_POS - current_POS) diff --git a/tests/cdb_reader/test_beam_strain_energy.py b/tests/cdb_reader/test_beam_strain_energy.py index da9f314..d23a2c2 100644 --- a/tests/cdb_reader/test_beam_strain_energy.py +++ b/tests/cdb_reader/test_beam_strain_energy.py @@ -18,15 +18,15 @@ _DATA = [ [101, 10, 1000, 0.34917488541448244], [102, 10, 1000, 1.2396276374589636], - [101, 10, 1001, 0.02940052665116039], + [101, 10, 1001, 0.029400526651160386], [102, 10, 1001, 0.0], - [101, 10, 1002, 0.3785754120656429], + [101, 10, 1002, 0.3785754120656428], [102, 10, 1002, 1.2396276374589636] ] _DATA_DIV = [ [101, 10, 1000, 0.3846355348887977], [102, 10, 1000, 0.8608525260131692], - [101, 10, 1001, 0.02940052665116039], + [101, 10, 1001, 0.029400526651160386], [102, 10, 1001, 0.0], [101, 10, 1002, 0.41403606153995803], [102, 10, 1002, 0.8608525260131692] @@ -73,7 +73,7 @@ def test_get_element_strain_energy(self) -> None: with self.subTest(msg="Existing entry"): self.assertEqual( self.cdb.beams.get_element_strain_energy(1002, 101), - 0.3785754120656429 + 0.3785754120656428 ) with self.subTest(msg="Non existing entry with default"): @@ -113,7 +113,7 @@ def test_get_strain_energy_after_clear_all(self) -> None: with self.subTest(msg="Check indexes management"): self.assertEqual( self.cdb.beams.get_element_strain_energy(1001, 101), - 0.02940052665116039 + 0.029400526651160386 ) @@ -197,5 +197,5 @@ def test_get_strain_energy_after_clear_all(self) -> None: with self.subTest(msg="Check indexes management"): self.assertEqual( self.cdb.beams.get_element_strain_energy(1001, 101), - 0.02940052665116039 + 0.029400526651160386 ) From 0a0161374e1ca456d60c6956829c50c4e84021e6 Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 12:24:51 +0200 Subject: [PATCH 13/16] assign cross-sectional values efficiently --- .../cdb_reader/_internals/beam.py | 26 ++++++++++++------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 292a7d6..afbd5da 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -60,7 +60,7 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: self.data.load() # Load properties for all unique PROP_END_1/2 - prop_nmb: list[int] = list( + prop_nmb = list( set(self.data.data()["PROP_END_1"].unique()) | set(self.data.data()["PROP_END_2"].unique()) ) @@ -81,15 +81,21 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: data["PROP_END_2"] = data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_2") for _ in elem_ids}).fillna(1).astype(int) data["LENGTH"] = data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) - for p in prop_nmb: - data.loc[data["PROP_END_1"] == p, "A"] = properties.get(p, "A") - data.loc[data["PROP_END_1"] == p, "IY"] = properties.get(p, "IY") - data.loc[data["PROP_END_1"] == p, "IZ"] = properties.get(p, "IZ") - data.loc[data["PROP_END_1"] == p, "EM"] = properties.get(p, "EM") - data.loc[data["PROP_END_2"] == p, "A"] = properties.get(p, "A") - data.loc[data["PROP_END_2"] == p, "IY"] = properties.get(p, "IY") - data.loc[data["PROP_END_2"] == p, "IZ"] = properties.get(p, "IZ") - data.loc[data["PROP_END_2"] == p, "EM"] = properties.get(p, "EM") + # Precompute property mappings + prop_mappings = { + "A": {_: properties.get(_, "A") for _ in prop_nmb}, + "IY": {_: properties.get(_, "IY") for _ in prop_nmb}, + "IZ": {_: properties.get(_, "IZ") for _ in prop_nmb}, + "EM": {_: properties.get(_, "EM") for _ in prop_nmb}, + } + + # Assign properties + for quantity in ["A", "IY", "IZ", "EM"]: + data[quantity] = data["PROP_END_1"].map(prop_mappings[quantity]) + mask_end_2 = data["POS_REL"] == 1 + data.loc[mask_end_2, quantity] = ( + data.loc[mask_end_2, "PROP_END_2"].map(prop_mappings[quantity]) + ) # Calculate strain energy data["U"] = data.eval( From d3c3571d79c22e7748a7d8b77c92d93fc8d6ff3c Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 12:32:35 +0200 Subject: [PATCH 14/16] avoid 2 x groupby by caching results --- src/py_sofistik_utils/cdb_reader/_internals/beam.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index afbd5da..501b9ae 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -102,14 +102,13 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: "(N ** 2 / A + MY ** 2 / IY + MZ ** 2 / IZ) / (2 * EM)" ) - # Calculate segment lengths (next_POS - current_POS) + # Calculate segment lengths data = data.sort_values(["LOAD_CASE", "ELEM_ID", "POS"]) - data["SEGMENT_LENGTH"] = data.groupby(["LOAD_CASE", "ELEM_ID"])["POS"].diff(-1).abs().fillna(0) - - # Calculate U for the next section (shifted U values) - data["U_NEXT"] = data.groupby(["LOAD_CASE", "ELEM_ID"])["U"].shift(-1).fillna(0) + grouped = data.groupby(["LOAD_CASE", "ELEM_ID"]) + data["SEGMENT_LENGTH"] = grouped["POS"].diff(-1).abs().fillna(0) # Calculate strain energy per segment: (U1 + U2) * L_seg / 2 + data["U_NEXT"] = grouped["U"].shift(-1).fillna(0) data["U_SEGMENT"] = (data["U"] + data["U_NEXT"]) * data["SEGMENT_LENGTH"] / 2 # Sum strain energy across all segments for each (LOAD_CASE, ELEM_ID) From b23016fd03f9c5230e04174516ea53acc42dc36d Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Tue, 5 May 2026 12:38:10 +0200 Subject: [PATCH 15/16] final clean-up --- .../cdb_reader/_internals/beam.py | 36 +++++++++++++------ 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 501b9ae..c0ba724 100644 --- a/src/py_sofistik_utils/cdb_reader/_internals/beam.py +++ b/src/py_sofistik_utils/cdb_reader/_internals/beam.py @@ -59,7 +59,7 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: if not self.data.is_loaded(): self.data.load() - # Load properties for all unique PROP_END_1/2 + # Load properties prop_nmb = list( set(self.data.data()["PROP_END_1"].unique()) | set(self.data.data()["PROP_END_2"].unique()) @@ -67,19 +67,31 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: properties = CrossSectionalData(self._dll) properties.load([int(_) for _ in prop_nmb]) - # Load results for all load cases + # Load results self.results.load(load_cases) - # Initialize data with zeros for required columns + # Initialize data data = self.results.data().reset_index(drop=True) - required_columns = ["LENGTH", "PROP_END_1", "PROP_END_2", "EM", "A", "IY", "IZ", "U"] + required_columns = [ + "LENGTH", "PROP_END_1", "PROP_END_2", "EM", "A", "IY", "IZ", "U" + ] data[required_columns] = 0.0 - # Assign PROP_END_1/2 and LENGTH using vectorized operations - elem_ids = data["ELEM_ID"].unique() - data["PROP_END_1"] = data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_1") for _ in elem_ids}).fillna(1).astype(int) - data["PROP_END_2"] = data["ELEM_ID"].map({_: self.data.get(_, "PROP_END_2") for _ in elem_ids}).fillna(1).astype(int) - data["LENGTH"] = data["ELEM_ID"].map({_: self.data.get(_) for _ in elem_ids}).fillna(1.0) + # Assign PROP_END_1/2 and LENGTH + ids = data["ELEM_ID"].unique() + data["PROP_END_1"] = ( + data["ELEM_ID"].map( + {_: self.data.get(_, "PROP_END_1") for _ in ids} + ).fillna(1).astype(int) + ) + data["PROP_END_2"] = ( + data["ELEM_ID"].map( + {_: self.data.get(_, "PROP_END_2") for _ in ids} + ).fillna(1).astype(int) + ) + data["LENGTH"] = ( + data["ELEM_ID"].map({_: self.data.get(_) for _ in ids}).fillna(1.0) + ) # Precompute property mappings prop_mappings = { @@ -109,7 +121,9 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: # Calculate strain energy per segment: (U1 + U2) * L_seg / 2 data["U_NEXT"] = grouped["U"].shift(-1).fillna(0) - data["U_SEGMENT"] = (data["U"] + data["U_NEXT"]) * data["SEGMENT_LENGTH"] / 2 + data["U_SEGMENT"] = ( + (data["U"] + data["U_NEXT"]) * data["SEGMENT_LENGTH"] / 2 + ) # Sum strain energy across all segments for each (LOAD_CASE, ELEM_ID) result = data.groupby(["LOAD_CASE", "ELEM_ID"]).agg({ @@ -117,7 +131,7 @@ def calculate_strain_energy(self, load_cases: int | list[int]) -> None: "GROUP": "first" }).rename(columns={"U_SEGMENT": "U"}).reset_index() - # Reorder columns to ['ELEM_ID', 'GROUP', 'LOAD_CASE', 'U'] + # Reorder columns result = result[["ELEM_ID", "GROUP", "LOAD_CASE", "U"]] # Set multi-index for fast lookups From 3cd2c2c527c75bb70689c0147887fac88baa61fc Mon Sep 17 00:00:00 2001 From: StudioWEngineers Date: Mon, 11 May 2026 14:16:15 +0200 Subject: [PATCH 16/16] add WIP docs --- docs/cdb/post/strain_energy.rst | 85 +++++++++++++++++++++++++++++++++ docs/cdb/post_processing.rst | 20 ++++++++ docs/index.rst | 6 +++ 3 files changed, 111 insertions(+) create mode 100644 docs/cdb/post/strain_energy.rst create mode 100644 docs/cdb/post_processing.rst diff --git a/docs/cdb/post/strain_energy.rst b/docs/cdb/post/strain_energy.rst new file mode 100644 index 0000000..a985eb1 --- /dev/null +++ b/docs/cdb/post/strain_energy.rst @@ -0,0 +1,85 @@ +Strain Energy +------------- + +The strain energy can be calculated from the internal forces obtained from +SOFiSTiK, together with the relevant cross-sectional properties and the +element length. + +Under the hyphoteses of: + +* small displacements and +* constant material and cross-sectional properties of the beam along its length, + +the general expression for the strain energy of an elastic beam is: + +.. math:: + :label: strain_energy_general + + U = \frac{1}{2} \int_{0}^{L} \left( + \frac{N^{2}\left(x\right)}{EA} + + \frac{M_{Y}^{2}\left(x\right)}{EI_{Y}} + + \frac{M_{Z}^{2}\left(x\right)}{EI_{Z}} + + \frac{k_{Y}V_{Y}^{2}\left(x\right)}{GA} + + \frac{k_{Z}V_{Z}^{2}\left(x\right)}{GA} + + \frac{T^{2}\left(x\right)}{GJ} + \right) dx + +where :math:`N` is the axial force, :math:`M_Y` and :math:`M_Z` are the +bending moments, :math:`V_Y` and :math:`V_Z` are the shear forces, :math:`T` is +the torsional moment, :math:`E` and :math:`G` are the normal and shear moduli, +:math:`A`, :math:`I_Y` and :math:`I_Z` are the area and the two second moments +of area of the cross-section. :math:`k_Y` and :math:`k_Z` are the two shear +factors. + +This expression is particularized and further developed for each type of finite +element in the next paraghraps. + +.. important:: + + The strain energy calculated by this module is approximate only. Further + details are provided below. + + Although efforts have been made to ensure that the relevant engineering + theory has been correctly implemented, it remains the user's + responsibility to verify and validate the results, particularly with + regard to whether the adopted simplifications are appropriate for the + intended use. + +Beam elements +""""""""""""" + +If the contributions of shear forces and torsion can be neglected, the +expression :eq:`strain_energy_general` can be simplified to: + +.. math:: + U = \frac{1}{2E} \int_{0}^{L} \left( + \frac{N^{2}\left(x\right)}{A} + + \frac{M_{Y}^{2}\left(x\right)}{I_{Y}} + + \frac{M_{Z}^{2}\left(x\right)}{I_{Z}} + \right) dx + +If we assume a generic linear variation for any of the internal force of the type: + +.. math:: + Q(x) = a + bx + +Its square integral over the beam length is: + +.. math:: + \int_{0}^{L} Q(x)^{2}dx = \frac{L}{3} \left( Q_{1}^{2} + Q_{1}Q_{2}+ Q_{2}^{2}\right) + +Where :math:`Q_1` and :math:`Q_2` are the values of :math:`Q` at beam end nodes + +Therefore the : + +.. math:: + + U = \frac{L}{6E} \left[ + \frac{N_{1}^{2} + N_{1}N_{2} + N_{2}^{2}}{A} + + + \frac{M_{Y1}^{2} + M_{Y1}M_{Y2} + M_{Y2}^{2}}{I_Y} + + + \frac{M_{Z1}^{2} + M_{Z1}M_{Z2} + M_{Z2}^{2}}{I_Z} + \right] + +This equation, which is the one implemented, is valid if the \ No newline at end of file diff --git a/docs/cdb/post_processing.rst b/docs/cdb/post_processing.rst new file mode 100644 index 0000000..d18e698 --- /dev/null +++ b/docs/cdb/post_processing.rst @@ -0,0 +1,20 @@ +.. _cdb_post_processing: + +Post-processing +=============== + +The CDBReader module provides functionality to calculate and/or manipulate +basic quantities extracted from the SOFiSTiK database, in addition to direct +access to raw data. + +Available functionality includes: + * calculation of strain energy + +The following subsections provide further details on this post-processing +capabilities. + +.. toctree:: + :maxdepth: 1 + :hidden: + + post/strain_energy diff --git a/docs/index.rst b/docs/index.rst index 0790f3e..fc8f666 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -53,6 +53,12 @@ for the terms and conditions of use. installation usage +.. toctree:: + :caption: CDB Reader + :hidden: + + cdb/post_processing + .. toctree:: :caption: Testing :maxdepth: 1