Update hinged component assembly, solver.

- Add back support for zerolength elements in solver.
- Add option to ignore elements based on their tag for a given analysis.

Author: ioannis-vm

src/osmg/analysis/load_case.py

@@ -18,14 +18,15 @@
 from osmg.analysis.supports import ElasticSupport, FixedSupport
 from osmg.core.common import EPSILON, THREE_DIMENSIONAL, TWO_DIMENSIONAL
 from osmg.core.osmg_collections import BeamColumnAssembly
-from osmg.model_objects.element import BeamColumnElement
+from osmg.model_objects.element import BeamColumnElement, Bar
 
 from tqdm import tqdm
 
 if TYPE_CHECKING:
     from osmg.analysis.common import PointLoad
     from osmg.core.model import Model, Model2D, Model3D
     from osmg.model_objects.node import Node
+    from osmg.core.osmg_collections import ComponentAssembly
 
 
 def ensure_minmax_level_exists_or_add(data: pd.DataFrame) -> pd.DataFrame:
@@ -241,7 +242,7 @@ def define_rigid_diaphragm(
     def calculate_basic_forces(  # noqa: C901
         self,
         recorder_name: str,
-        element_lengths: dict[int, float],
+        components: dict[int, ComponentAssembly],
         *,
         ndm: int,
         num_stations: int = 12,
@@ -300,9 +301,25 @@ def calculate_basic_forces(  # noqa: C901
             raise TypeError(msg)
 
         if isinstance(self, HasLoads):
-            udls = self.load_registry.element_udl
+            udls_global = self.load_registry.element_udl
         else:
-            udls = {}
+            udls_global = {}
+
+        udls_local = {}
+        elements = {}
+        for component_uid, component in components.items():
+            if not isinstance(component, BeamColumnAssembly):
+                continue
+            for element in component.elements.values():
+                if isinstance(element, (BeamColumnElement, Bar)):
+                    elements[element.uid] = element
+            elements.update(component.elements)
+            component_global_udl = udls_global.get(component_uid)
+            if component_global_udl:
+                component_local_udls = component.calculate_element_udl(
+                    component_global_udl
+                )
+                udls_local.update(component_local_udls)
 
         data = recorder.get_data()
 
@@ -327,11 +344,17 @@ def calculate_basic_forces(  # noqa: C901
         #     dof: df.xs(1.0, level='station', axis=1) for dof, df in dof_data.items()
         # }
 
-        missing_elements = [e for e in recorder.elements if e not in element_lengths]
+        missing_elements = [e for e in recorder.elements if e not in elements]
         if missing_elements:
-            msg = f'Missing lengths for elements: {missing_elements}'
+            msg = f'Missing elements: {missing_elements}'
             raise ValueError(msg)
 
+        element_lengths: dict[int, float] = {
+            element.uid: element.clear_length()
+            for element in elements.values()
+            if isinstance(element, BeamColumnElement)
+        }
+
         locations = np.array(
             [
                 np.linspace(0.00, element_lengths[element], num=num_stations)
@@ -360,7 +383,9 @@ def calculate_basic_forces(  # noqa: C901
         # data_j_axial = dof_data_j[1]
         w_data_axial = np.array(
             [
-                udls[element_uid][0] if udls.get(element_uid) is not None else 0.00
+                udls_local[element_uid][0]
+                if udls_local.get(element_uid) is not None
+                else 0.00
                 for element_uid in recorder.elements
             ]
         )[np.newaxis, :, np.newaxis]
@@ -385,7 +410,9 @@ def calculate_basic_forces(  # noqa: C901
         # data_j_shear_y = dof_data_j[2]
         w_data_shear_y = np.array(
             [
-                udls[element_uid][1] if udls.get(element_uid) is not None else 0.00
+                udls_local[element_uid][1]
+                if udls_local.get(element_uid) is not None
+                else 0.00
                 for element_uid in recorder.elements
             ]
         )[np.newaxis, :, np.newaxis]
@@ -418,8 +445,8 @@ def calculate_basic_forces(  # noqa: C901
             # data_j_shear_z = dof_data_j[3]
             w_data_shear_z = np.array(
                 [
-                    udls[element_uid][2]
-                    if udls.get(element_uid) is not None
+                    udls_local[element_uid][2]
+                    if udls_local.get(element_uid) is not None
                     else 0.00
                     for element_uid in recorder.elements
                 ]
@@ -680,18 +707,16 @@ def self_weight(self, case_name: str, scaling_factor: float = 1.0) -> None:
           case_name: Name of the load case to be created.
           scaling_factor: Self-weight scaling factor to use.
         """
-        # get all beamcolumn elements
-        elements = [
-            element
+        # get all beamcolumn assemblies
+        components = [
+            component
             for component in self.model.components.values()
             if isinstance(component, BeamColumnAssembly)
-            for element in component.elements.values()
-            if isinstance(element, BeamColumnElement)
         ]
-        for element in elements:
-            weight_per_length = element.section.sec_w * scaling_factor
+        for component in components:
+            weight_per_length = component.get_section().sec_w * scaling_factor
             udl = UDL((0.00, 0.00, -weight_per_length))
-            self.dead[case_name].load_registry.element_udl[element.uid] = udl
+            self.dead[case_name].load_registry.element_udl[component.uid] = udl
 
     def self_mass(
         self,
@@ -712,30 +737,29 @@ def self_mass(
                              a scaling coefficient.
           g_constant: Factor to convert loads to mass.
         """
-        # get all beamcolumn elements
-        elements = {
-            element.uid: element
+        # get all BeamColumn assemblies
+        components = {
+            component.uid: component
             for component in self.model.components.values()
             if isinstance(component, BeamColumnAssembly)
-            for element in component.elements.values()
-            if isinstance(element, BeamColumnElement)
         }
         for source_load_case, factor in source_load_cases:
             # Convert UDL to mass
             for uid, udl in source_load_case.load_registry.element_udl.items():
-                element = elements[uid]
-                length = element.clear_length()
+                component = components[uid]
+                length = component.clear_length()
                 weight = np.abs(udl[-1] * length)
-                mass = weight * factor / g_constant / 2.0
+                num_nodes = len(component.internal_nodes)
+                mass = weight * factor / g_constant / num_nodes
                 # TODO(JVM): separate cases for other ndm/ndf
                 # configurations.
                 point_mass = PointMass((mass, mass, mass, 0.00, 0.00, 0.00))
-                for node in element.nodes:
-                    if node.uid not in target_load_case.mass_registry:
-                        target_load_case.mass_registry[node.uid] = point_mass
+                for node_uid in component.internal_nodes:
+                    if node_uid not in target_load_case.mass_registry:
+                        target_load_case.mass_registry[node_uid] = point_mass
                     else:
-                        existing_mass = target_load_case.mass_registry[node.uid]
-                        target_load_case.mass_registry[node.uid] = PointMass(
+                        existing_mass = target_load_case.mass_registry[node_uid]
+                        target_load_case.mass_registry[node_uid] = PointMass(
                             (*(e + p for e, p in zip(existing_mass, point_mass)),)
                         )
 

src/osmg/analysis/solver.py

@@ -22,6 +22,7 @@
     DispBeamColumn,
     ElasticBeamColumn,
     TwoNodeLink,
+    ZeroLength,
 )
 
 try:
@@ -49,6 +50,7 @@ class AnalysisSettings:
     constraints: tuple[str] = ('Transformation',)
     numberer: Literal['RCM'] = 'RCM'
     system: Literal['Umfpack'] = 'Umfpack'
+    ignore_by_tag: set[str] = field(default_factory=set)
 
 
 @dataclass()
@@ -152,10 +154,9 @@ def opensees_define_material(self, material: UniaxialMaterial) -> None:
             ops.uniaxialMaterial(*material.ops_args())
             self._defined_materials.append(material.uid)
 
-    @staticmethod
-    def opensees_define_nodes(model: Model) -> None:
+    def opensees_define_nodes(self, model: Model) -> None:
         """Define the nodes of the model in OpenSees."""
-        for uid, node in model.get_all_nodes().items():
+        for uid, node in model.get_all_nodes(self.settings.ignore_by_tag).items():
             ops.node(uid, *node.coordinates)
 
     def opensees_define_elements(
@@ -166,6 +167,7 @@ def opensees_define_elements(
         elastic_beamcolumn_elements: list[ElasticBeamColumn] = []
         bar_elements: list[Bar] = []
         two_node_link_elements: list[TwoNodeLink] = []
+        zerolength_elements: list[ZeroLength] = []
         unsupported_element_types: list[str] = []
 
         # Note: Materials are defined on an as-needed basis.  We keep
@@ -176,6 +178,8 @@ def opensees_define_elements(
 
         components = model.components.values()
         for component in components:
+            if component.tags & self.settings.ignore_by_tag:
+                continue
             elements = component.elements
             for element in elements.values():
                 if isinstance(element, ElasticBeamColumn):
@@ -184,6 +188,8 @@ def opensees_define_elements(
                     bar_elements.append(element)
                 elif isinstance(element, TwoNodeLink):
                     two_node_link_elements.append(element)
+                elif isinstance(element, ZeroLength):
+                    zerolength_elements.append(element)
                 else:
                     unsupported_element_types.append(element.__class__.__name__)
 
@@ -197,6 +203,7 @@ def opensees_define_elements(
         )
         self.opensees_define_bar_elements(bar_elements)
         self.opensees_define_two_node_link_elements(two_node_link_elements)
+        self.opensees_define_zerolength_elements(zerolength_elements)
 
         # clear defined materials
         self._defined_materials = []
@@ -239,6 +246,15 @@ def opensees_define_two_node_link_elements(
                 self.opensees_define_material(material)
             ops.element(*element.ops_args())
 
+    def opensees_define_zerolength_elements(
+        self, elements: list[ZeroLength]
+    ) -> None:
+        """Define Zerolength elements."""
+        for element in elements:
+            for material in element.materials:
+                self.opensees_define_material(material)
+            ops.element(*element.ops_args())
+
     def opensees_define_node_restraints(
         self, model: Model, load_case: LoadCase
     ) -> None:
@@ -299,8 +315,8 @@ def opensees_define_model(self, model: Model, load_case: LoadCase) -> None:
             self.define_default_recorders(model)
         self.opensees_define_recorders()
 
-    @staticmethod
     def opensees_define_loads(
+        self,
         model: Model,
         load_case: LoadCase,
         amplification_factor: float = 1.00,
@@ -326,6 +342,8 @@ def opensees_define_loads(
         # UDL on components
         for component_uid, global_udl in load_case.load_registry.element_udl.items():  # type: ignore
             component = model.components[component_uid]
+            if component.tags & self.settings.ignore_by_tag:
+                continue
             assert isinstance(component, BeamColumnAssembly)
             local_udls = component.calculate_element_udl(global_udl)
             for beamcolumn_element_uid, local_udl in local_udls.items():
@@ -389,6 +407,8 @@ def define_default_recorders(self, model: Model) -> None:
         applicable_elements = []
         components = model.components.values()
         for component in components:
+            if component.tags & self.settings.ignore_by_tag:
+                continue
             elements = component.elements
             for element in elements.values():
                 if isinstance(element, (ElasticBeamColumn, DispBeamColumn, Bar)):
@@ -612,15 +632,9 @@ def run(self, model: Model, load_case: LoadCase) -> None:  # noqa: C901  # type:
             # in the dict.
 
             # Verifying displacements are correct.
-            for some_node in mode_eigenvectors.index.get_level_values(
-                'node'
-            ).unique():
-                for check_dof in range(1, ndf + 1):
-                    if (node, check_dof) in mode_eigenvectors:
-                        assert np.allclose(
-                            np.array(ops.nodeDisp(some_node))[check_dof - 1],
-                            mode_eigenvectors[some_node, check_dof],
-                        )
+            for key, value in mode_eigenvectors.items():
+                node_uid, dof = key
+                assert np.allclose(ops.nodeDisp(node_uid, dof), value)
 
             self.log('   Wiping OpenSees domain.')
             # Doing this before reading the basic force recorder data

src/osmg/core/model.py

@@ -69,10 +69,15 @@ def reference_length(self) -> float:
         p_min, p_max = self.bounding_box(padding=0.00)
         return float(np.max(p_max - p_min))
 
-    def get_all_nodes(self) -> dict[int, Node]:
+    def get_all_nodes(
+        self, ignore_by_tag: set[str] | None = None
+    ) -> dict[int, Node]:
         """
         Get all nodes in the model.
 
+        Params:
+          ignore_by_tag: Set of tags of components to ignore.
+
         Returns:
           A dictionary with the nodes. Keys are their UIDs.
         """
@@ -82,6 +87,8 @@ def get_all_nodes(self) -> dict[int, Node]:
         internal_nodes: dict[int, Node] = {}
         components = self.components.values()
         for component in components:
+            if ignore_by_tag and component.tags & ignore_by_tag:
+                continue
             internal_nodes.update(component.internal_nodes)
 
         all_nodes: dict[int, Node] = {}