feat: add option to calculate the complete mn domain

docs/usage/sections/index.md

@@ -54,6 +54,11 @@ Furthermore, we have the following methods:
 
 See the {class}`GenericSectionCalculator <structuralcodes.sections.GenericSectionCalculator>` for a complete list.
 
+(usage-sections-complete-nm-note)=
+:::{note}
+Notice that a call to {py:meth}`.calculate_nm_interaction_domain() <structuralcodes.sections.GenericSectionCalculator.calculate_nm_interaction_domain>` returns the interaction domain for a negative bending moment, i.e. a bending moment that gives compression at the top of the cross section according to the [sign convention](#theory-sign-convention). To obtain the interaction domain for the positive bending moment, the neutral axis for the calculation should be rotated an angle {math}`\theta = \pi`, i.e. the method should be called with the keyword argument `theta = np.pi` or `theta = math.pi`. Alternatively, to return the complete domain, the method could be called with the keyword argument `complete_domain = True`.
+:::
+
 (code-usage-generic-section)=
 ::::{dropdown-syntax}
 :::{literalinclude} ../../_example_code/usage_create_surface_geometries_with_reinforcement.py

structuralcodes/sections/_generic.py

@@ -1003,6 +1003,7 @@ def calculate_nm_interaction_domain(
         type_4: t.Literal['linear', 'geometric', 'quadratic'] = 'linear',
         type_5: t.Literal['linear', 'geometric', 'quadratic'] = 'linear',
         type_6: t.Literal['linear', 'geometric', 'quadratic'] = 'linear',
+        complete_domain: bool = False,
     ) -> s_res.NMInteractionDomain:
         """Calculate the NM interaction domain.
 
@@ -1037,6 +1038,9 @@ def calculate_nm_interaction_domain(
                 type_1 for options.
             type_6 (str): Type of spacing for field 6 (default = 'linear'). See
                 type_1 for options.
+            complete_domain (bool): Flag to indicate if only the part of the
+                domain related to the negative moment is returned (False), or
+                if the positive part of the domain is also included (True).
 
         Returns:
             NMInteractionDomain: The calculation results.
@@ -1053,7 +1057,7 @@ def calculate_nm_interaction_domain(
                 )
             )
 
-        # Get ultimate strain profiles for theta angle
+        # Get ultimate strain profiles for theta = theta
         strains, field_num = self._compute_ultimate_strain_profiles(
             theta=theta,
             num_1=num_1,
@@ -1070,6 +1074,30 @@ def calculate_nm_interaction_domain(
             type_6=type_6,
         )
 
+        if complete_domain:
+            # Add strain profiles for theta = (theta + pi)
+            additional_strains, additional_field_num = (
+                self._compute_ultimate_strain_profiles(
+                    theta=theta + np.pi,
+                    num_1=num_1,
+                    num_2=num_2,
+                    num_3=num_3,
+                    num_4=num_4,
+                    num_5=num_5,
+                    num_6=num_6,
+                    type_1=type_1,
+                    type_2=type_2,
+                    type_3=type_3,
+                    type_4=type_4,
+                    type_5=type_5,
+                    type_6=type_6,
+                )
+            )
+            strains = np.concatenate((strains, additional_strains[-2:0:-1]))
+            field_num = np.concatenate(
+                (field_num, additional_field_num[-2:0:-1])
+            )
+
         # integrate all strain profiles
         forces = np.zeros_like(strains)
         for i, strain in enumerate(strains):

tests/test_sections/test_generic_section.py

@@ -1658,3 +1658,79 @@ def test_issue_cracked_properties():
     assert cracked_properties_before.isclose(
         cracked_properties_after, rtol=1e-3, atol=1e-6
     )
+
+
+def test_mn_full_domain():
+    """Test calculating the full MN interaction domain."""
+    # Set parameters
+    width = 250
+    height = 500
+    diameter_reinf = 25
+    cover = 50
+
+    fck = 45
+    fyk = 500
+    Es = 200e3
+    epsuk = 6e-2
+
+    # Create materials
+    concrete = ConcreteEC2_2004(fck=fck)
+    reinforcement = ReinforcementEC2_2004(fyk=fyk, Es=Es, ftk=fyk, epsuk=epsuk)
+
+    # Create geometry
+    z_reinforcement = height / 2 - cover - diameter_reinf / 2
+    y_reinforcement = width / 2 - cover - diameter_reinf / 2
+    geometry = RectangularGeometry(
+        width=width, height=height, material=concrete
+    )
+
+    for z in (-z_reinforcement, z_reinforcement):
+        geometry = add_reinforcement_line(
+            geometry,
+            (z, -y_reinforcement),
+            (z, y_reinforcement),
+            diameter_reinf,
+            reinforcement,
+            n=2,
+        )
+
+    # Create section
+    section = GenericSection(geometry, integrator='fiber')
+
+    # Calculate interaction domain for theta = 0
+    interaction_domain_0 = (
+        section.section_calculator.calculate_nm_interaction_domain(theta=0)
+    )
+
+    # Calculate interaction domain for theta = pi
+    interaction_domain_180 = (
+        section.section_calculator.calculate_nm_interaction_domain(theta=np.pi)
+    )
+
+    # Calculate the full interaction domain
+    interaction_domain_full = (
+        section.section_calculator.calculate_nm_interaction_domain(
+            complete_domain=True
+        )
+    )
+
+    # Combine theta = 0 and theta = 180 to obtain the full domain
+    interaction_domain_full_combined_n = [
+        *interaction_domain_0.n,
+        *interaction_domain_180.n[-2:0:-1],
+    ]
+    interaction_domain_full_combined_my = [
+        *interaction_domain_0.m_y,
+        *interaction_domain_180.m_y[-2:0:-1],
+    ]
+
+    assert (
+        len(interaction_domain_full.n)
+        == len(interaction_domain_0.n) + len(interaction_domain_180.n) - 2
+    )
+    assert np.allclose(
+        interaction_domain_full.n, interaction_domain_full_combined_n
+    )
+    assert np.allclose(
+        interaction_domain_full.m_y, interaction_domain_full_combined_my
+    )