moment_interaction_diagram() crashes Python

[jindrichpetrasek-alt]

Hi @robbievanleeuwen,
I'm having issues calculating the M-N interaction diagram for a concrete section with a stress-strain profile other than rectangular block. For some unknown reason, the following code crashes Python without generating any error. As far as I can tell, the problem might be related to meshing or cytriangle? Or is it just a problem between keyboard and chair and there is something wrong with my input?

import numpy as np
from sectionproperties.pre.library.concrete_sections import concrete_rectangular_section
from sectionproperties.pre. geometry import Geometry

import concreteproperties.stress_strain_profile as ssp
from concreteproperties.material import Concrete, SteelBar
from concreteproperties.concrete_section import ConcreteSection
from concreteproperties.results import MomentInteractionResults
from concreteproperties.post import si_kn_m

concrete = Concrete(
    name="C30/37",
    density=2.5e-6,
    stress_strain_profile=ssp.ConcreteLinearNoTension(elastic_modulus=3.284000e+04,
                                                         ultimate_strain=3.500000e-03,
                                                         compressive_strength=38.00),
    ultimate_stress_strain_profile=ssp.BilinearStressStrain(compressive_strength=30/1.5,
                                                compressive_strain=1.7500000e-03,
                                                ultimate_strain=3.500000e-03 ),
    flexural_tensile_strength=2.9,
    colour="lightgrey")

steel = SteelBar(
    name="B500B",
    density=7.85e-6,
    stress_strain_profile=ssp.SteelElasticPlastic( 
        yield_strength=500/1.15,
        elastic_modulus=200e3,
        fracture_strain=0.05),
    colour="grey")

geom = concrete_rectangular_section(
    d=900,
    b=300,
    dia_top=16,
    area_top=np.pi*8*8,
    n_top=3,
    c_top=42,
    dia_bot=16,
    area_bot=np.pi*8*8,
    n_bot=3,
    c_bot=42,
    conc_mat=concrete,
    steel_mat=steel)

conc_sec = ConcreteSection(geom)
conc_sec.plot_section()

res_ultM = conc_sec.ultimate_bending_capacity(n=3000*1e3)
print(f"MRd (N = -3000 kN): {res_ultM.m_xy / 1e6:.2f} kN.m")

res_intNM = []
res_intNM.append(conc_sec.moment_interaction_diagram(theta=0, progress_bar=True, n_points=12))
res_intNM.append(conc_sec.moment_interaction_diagram(theta=np.pi, progress_bar=True, n_points=12))


MomentInteractionResults.plot_multiple_diagrams(
    moment_interaction_results=res_intNM,
    labels=["Positive", "Negative"],
    fmt="-",
    units=si_kn_m)

Thank you

[github-actions]

Thanks for opening your first issue in concreteproperties 🙌 Pull requests are always welcome 😉

[ocabitza]

Hi @jindrichpetrasek-alt,
I've managed to make it work by simply modifying the theta angle of the negative plot: if you set theta = -np.pi instead of np.pi, it works.

I don't know why it crashes with theta = np.pi. As you already suggested, the bug seems to be related to cytriangle and meshing (the program raises an error if you set theta = 0.99*np.pi, as you have probably already noticed).