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/cdb/test_setup.rst b/docs/cdb/test_setup.rst index c8c8982..f4fb27d 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..a338e63 --- /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 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 diff --git a/src/py_sofistik_utils/cdb_reader/_internals/beam.py b/src/py_sofistik_utils/cdb_reader/_internals/beam.py index 430d0a7..c0ba724 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 . cross_section_data import CrossSectionalData from . sofistik_dll import SofDll @@ -31,5 +32,182 @@ 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 + prop_nmb = list( + set(self.data.data()["PROP_END_1"].unique()) | + set(self.data.data()["PROP_END_2"].unique()) + ) + properties = CrossSectionalData(self._dll) + properties.load([int(_) for _ in prop_nmb]) + + # Load results + self.results.load(load_cases) + + # Initialize data + data = self.results.data().reset_index(drop=True) + 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 + 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 = { + "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( + "(N ** 2 / A + MY ** 2 / IY + MZ ** 2 / IZ) / (2 * EM)" + ) + + # Calculate segment lengths + data = data.sort_values(["LOAD_CASE", "ELEM_ID", "POS"]) + 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) + result = data.groupby(["LOAD_CASE", "ELEM_ID"]).agg({ + "U_SEGMENT": "sum", + "GROUP": "first" + }).rename(columns={"U_SEGMENT": "U"}).reset_index() + + # Reorder columns + 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) 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 e8af549..d62a313 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 get( self, @@ -147,6 +149,10 @@ 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 get_data(self, deep: bool = True) -> DataFrame: """Return the :class:`pandas.DataFrame` containing the loaded key ``100/00``. @@ -276,3 +282,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 diff --git a/tests/cdb_reader/test_beam_data.py b/tests/cdb_reader/test_beam_data.py index 3d04b9f..3839796 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()) 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..d23a2c2 --- /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.029400526651160386], + [102, 10, 1001, 0.0], + [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.029400526651160386], + [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.3785754120656428 + ) + + 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.029400526651160386 + ) + + +@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.029400526651160386 + )