from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from typing import TYPE_CHECKING, List, Optional, Tuple, Union
import matplotlib.cm as cm
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
from matplotlib.collections import PatchCollection
from matplotlib.colors import CenteredNorm # type: ignore
from rich.console import Console
from rich.table import Table
from scipy.interpolate import interp1d
from sectionproperties.pre.geometry import CompoundGeometry
from shapely import Point, Polygon
from concreteproperties.post import plotting_context
if TYPE_CHECKING:
import matplotlib
from concreteproperties.analysis_section import AnalysisSection
from concreteproperties.concrete_section import ConcreteSection
from concreteproperties.pre import CPGeom
[docs]@dataclass
class GrossProperties:
"""Class for storing gross concrete section properties.
All properties with an `e_` preceding the property are multiplied by the elastic
modulus. In order to obtain transformed properties, call the
:meth:`~concreteproperties.concrete_section.ConcreteSection.get_transformed_gross_properties`
method.
"""
# section areas
total_area: float = 0
concrete_area: float = 0
reinf_meshed_area: float = 0
reinf_lumped_area: float = 0
strand_area: float = 0
e_a: float = 0
# section mass
mass: float = 0
# section perimeter
perimeter: float = 0
# first moments of area
e_qx: float = 0
e_qy: float = 0
qx_gross: float = 0
qy_gross: float = 0
# centroids
cx: float = 0
cy: float = 0
cx_gross: float = 0
cy_gross: float = 0
# second moments of area
e_ixx_g: float = 0
e_iyy_g: float = 0
e_ixy_g: float = 0
e_ixx_c: float = 0
e_iyy_c: float = 0
e_ixy_c: float = 0
e_i11: float = 0
e_i22: float = 0
# principal axis angle
phi: float = 0
# section moduli
e_zxx_plus: float = 0
e_zxx_minus: float = 0
e_zyy_plus: float = 0
e_zyy_minus: float = 0
e_z11_plus: float = 0
e_z11_minus: float = 0
e_z22_plus: float = 0
e_z22_minus: float = 0
# other properties
conc_ultimate_strain: float = 0
n_prestress: float = 0
m_prestress: float = 0
[docs] def print_results(
self,
fmt: str = "8.6e",
):
"""Prints the gross concrete section properties to the terminal.
:param fmt: Number format
"""
table = Table(title="Gross Concrete Section Properties")
table.add_column("Property", justify="left", style="cyan", no_wrap=True)
table.add_column("Value", justify="right", style="green")
table.add_row("Total Area", "{:>{fmt}}".format(self.total_area, fmt=fmt))
table.add_row("Concrete Area", "{:>{fmt}}".format(self.concrete_area, fmt=fmt))
if self.reinf_meshed_area:
table.add_row(
"Meshed Reinforcement Area",
"{:>{fmt}}".format(self.reinf_meshed_area, fmt=fmt),
)
table.add_row(
"Lumped Reinforcement Area",
"{:>{fmt}}".format(self.reinf_lumped_area, fmt=fmt),
)
if self.strand_area:
table.add_row("Strand Area", "{:>{fmt}}".format(self.strand_area, fmt=fmt))
table.add_row("Axial Rigidity (EA)", "{:>{fmt}}".format(self.e_a, fmt=fmt))
table.add_row("Mass (per unit length)", "{:>{fmt}}".format(self.mass, fmt=fmt))
table.add_row("Perimeter", "{:>{fmt}}".format(self.perimeter, fmt=fmt))
table.add_row("E.Qx", "{:>{fmt}}".format(self.e_qx, fmt=fmt))
table.add_row("E.Qy", "{:>{fmt}}".format(self.e_qy, fmt=fmt))
table.add_row("x-Centroid", "{:>{fmt}}".format(self.cx, fmt=fmt))
table.add_row("y-Centroid", "{:>{fmt}}".format(self.cy, fmt=fmt))
table.add_row("x-Centroid (Gross)", "{:>{fmt}}".format(self.cx_gross, fmt=fmt))
table.add_row("y-Centroid (Gross)", "{:>{fmt}}".format(self.cy_gross, fmt=fmt))
table.add_row("E.Ixx_g", "{:>{fmt}}".format(self.e_ixx_g, fmt=fmt))
table.add_row("E.Iyy_g", "{:>{fmt}}".format(self.e_iyy_g, fmt=fmt))
table.add_row("E.Ixy_g", "{:>{fmt}}".format(self.e_ixy_g, fmt=fmt))
table.add_row("E.Ixx_c", "{:>{fmt}}".format(self.e_ixx_c, fmt=fmt))
table.add_row("E.Iyy_c", "{:>{fmt}}".format(self.e_iyy_c, fmt=fmt))
table.add_row("E.Ixy_c", "{:>{fmt}}".format(self.e_ixy_c, fmt=fmt))
table.add_row("E.I11", "{:>{fmt}}".format(self.e_i11, fmt=fmt))
table.add_row("E.I22", "{:>{fmt}}".format(self.e_i22, fmt=fmt))
table.add_row("Principal Axis Angle", "{:>{fmt}}".format(self.phi, fmt=fmt))
table.add_row("E.Zxx+", "{:>{fmt}}".format(self.e_zxx_plus, fmt=fmt))
table.add_row("E.Zxx-", "{:>{fmt}}".format(self.e_zxx_minus, fmt=fmt))
table.add_row("E.Zyy+", "{:>{fmt}}".format(self.e_zyy_plus, fmt=fmt))
table.add_row("E.Zyy-", "{:>{fmt}}".format(self.e_zyy_minus, fmt=fmt))
table.add_row("E.Z11+", "{:>{fmt}}".format(self.e_z11_plus, fmt=fmt))
table.add_row("E.Z11-", "{:>{fmt}}".format(self.e_z11_minus, fmt=fmt))
table.add_row("E.Z22+", "{:>{fmt}}".format(self.e_z22_plus, fmt=fmt))
table.add_row("E.Z22-", "{:>{fmt}}".format(self.e_z22_minus, fmt=fmt))
table.add_row(
"Ultimate Concrete Strain",
"{:>{fmt}}".format(self.conc_ultimate_strain, fmt=fmt),
)
# add prestressed results if they exist
if self.n_prestress:
table.add_row("n_prestress", "{:>{fmt}}".format(self.n_prestress, fmt=fmt))
table.add_row("m_prestress", "{:>{fmt}}".format(self.m_prestress, fmt=fmt))
console = Console()
console.print(table)
[docs]@dataclass
class CrackedResults:
r"""Class for storing cracked concrete section properties.
All properties with an `e_` preceding the property are multiplied by the elastic
modulus. In order to obtain transformed properties, call the
:meth:`~concreteproperties.results.CrackedResults.calculate_transformed_properties`
method.
:param theta: Angle (in radians) the neutral axis makes with the horizontal
axis (:math:`-\pi \leq \theta \leq \pi`)
"""
theta: float
n: float = 0
m: float = 0
m_cr: Union[float, Tuple[float, float]] = 0
d_nc: float = 0
cracked_geometries: List[CPGeom] = field(default_factory=list, repr=False)
e_a_cr: float = 0
e_qx_cr: float = 0
e_qy_cr: float = 0
cx: float = 0
cy: float = 0
e_ixx_g_cr: float = 0
e_iyy_g_cr: float = 0
e_ixy_g_cr: float = 0
e_ixx_c_cr: float = 0
e_iyy_c_cr: float = 0
e_ixy_c_cr: float = 0
e_iuu_cr: float = 0
e_i11_cr: float = 0
e_i22_cr: float = 0
phi_cr: float = 0
# transformed properties
elastic_modulus_ref: Optional[float] = None
a_cr: Optional[float] = None
qx_cr: Optional[float] = None
qy_cr: Optional[float] = None
ixx_g_cr: Optional[float] = None
iyy_g_cr: Optional[float] = None
ixy_g_cr: Optional[float] = None
ixx_c_cr: Optional[float] = None
iyy_c_cr: Optional[float] = None
ixy_c_cr: Optional[float] = None
iuu_cr: Optional[float] = None
i11_cr: Optional[float] = None
i22_cr: Optional[float] = None
[docs] def reset_results(self) -> None:
"""Resets the analysis results."""
self.e_a_cr = 0
self.e_qx_cr = 0
self.e_qy_cr = 0
self.cx = 0
self.cy = 0
self.e_ixx_g_cr = 0
self.e_iyy_g_cr = 0
self.e_ixy_g_cr = 0
self.e_ixx_c_cr = 0
self.e_iyy_c_cr = 0
self.e_ixy_c_cr = 0
self.e_iuu_cr = 0
self.e_i11_cr = 0
self.e_i22_cr = 0
self.phi_cr = 0
[docs] def plot_cracked_geometries(
self,
title: str = "Cracked Geometries",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots the geometries that remain (are in compression or are reinforcement)
after a cracked analysis.
:param title: Plot title
:param kwargs: Passed to
:meth:`~sectionproperties.pre.geometry.CompoundGeometry.plot_geometry`
:return: Matplotlib axes object
"""
return CompoundGeometry(
[geom.to_sp_geom() for geom in self.cracked_geometries]
).plot_geometry(title=title, **kwargs)
[docs] def print_results(
self,
fmt: str = "8.6e",
):
"""Prints the cracked concrete section properties to the terminal.
:param fmt: Number format
"""
table = Table(title="Cracked Concrete Section Properties")
table.add_column("Property", justify="left", style="cyan", no_wrap=True)
table.add_column("Value", justify="right", style="green")
table.add_row("theta", "{:>{fmt}}".format(self.theta, fmt=fmt))
table.add_row("n", "{:>{fmt}}".format(self.n, fmt=fmt))
table.add_row("m", "{:>{fmt}}".format(self.m, fmt=fmt))
if self.elastic_modulus_ref:
table.add_row(
"E_ref", "{:>{fmt}}".format(self.elastic_modulus_ref, fmt=fmt)
)
if isinstance(self.m_cr, tuple):
table.add_row("m_cr_pos", "{:>{fmt}}".format(self.m_cr[0], fmt=fmt))
table.add_row("m_cr_neg", "{:>{fmt}}".format(self.m_cr[1], fmt=fmt))
else:
table.add_row("m_cr", "{:>{fmt}}".format(self.m_cr, fmt=fmt))
table.add_row("d_nc", "{:>{fmt}}".format(self.d_nc, fmt=fmt))
if self.a_cr:
table.add_row("A_cr", "{:>{fmt}}".format(self.a_cr, fmt=fmt))
table.add_row("E.A_cr", "{:>{fmt}}".format(self.e_a_cr, fmt=fmt))
if self.qx_cr:
table.add_row("Qx_cr", "{:>{fmt}}".format(self.qx_cr, fmt=fmt))
table.add_row("Qy_cr", "{:>{fmt}}".format(self.qy_cr, fmt=fmt))
table.add_row("E.Qx_cr", "{:>{fmt}}".format(self.e_qx_cr, fmt=fmt))
table.add_row("E.Qy_cr", "{:>{fmt}}".format(self.e_qy_cr, fmt=fmt))
table.add_row("x-Centroid", "{:>{fmt}}".format(self.cx, fmt=fmt))
table.add_row("y-Centroid", "{:>{fmt}}".format(self.cy, fmt=fmt))
if self.ixx_g_cr:
table.add_row("Ixx_g_cr", "{:>{fmt}}".format(self.ixx_g_cr, fmt=fmt))
table.add_row("Iyy_g_cr", "{:>{fmt}}".format(self.iyy_g_cr, fmt=fmt))
table.add_row("Ixy_g_cr", "{:>{fmt}}".format(self.ixy_g_cr, fmt=fmt))
table.add_row("Ixx_c_cr", "{:>{fmt}}".format(self.ixx_c_cr, fmt=fmt))
table.add_row("Iyy_c_cr", "{:>{fmt}}".format(self.iyy_c_cr, fmt=fmt))
table.add_row("Ixy_c_cr", "{:>{fmt}}".format(self.ixy_c_cr, fmt=fmt))
table.add_row("Iuu_cr", "{:>{fmt}}".format(self.iuu_cr, fmt=fmt))
table.add_row("I11_cr", "{:>{fmt}}".format(self.i11_cr, fmt=fmt))
table.add_row("I22_cr", "{:>{fmt}}".format(self.i22_cr, fmt=fmt))
table.add_row("E.Ixx_g_cr", "{:>{fmt}}".format(self.e_ixx_g_cr, fmt=fmt))
table.add_row("E.Iyy_g_cr", "{:>{fmt}}".format(self.e_iyy_g_cr, fmt=fmt))
table.add_row("E.Ixy_g_cr", "{:>{fmt}}".format(self.e_ixy_g_cr, fmt=fmt))
table.add_row("E.Ixx_c_cr", "{:>{fmt}}".format(self.e_ixx_c_cr, fmt=fmt))
table.add_row("E.Iyy_c_cr", "{:>{fmt}}".format(self.e_iyy_c_cr, fmt=fmt))
table.add_row("E.Ixy_c_cr", "{:>{fmt}}".format(self.e_ixy_c_cr, fmt=fmt))
table.add_row("E.Iuu_cr", "{:>{fmt}}".format(self.e_iuu_cr, fmt=fmt))
table.add_row("E.I11_cr", "{:>{fmt}}".format(self.e_i11_cr, fmt=fmt))
table.add_row("E.I22_cr", "{:>{fmt}}".format(self.e_i22_cr, fmt=fmt))
table.add_row("phi_cr", "{:>{fmt}}".format(self.phi_cr, fmt=fmt))
console = Console()
console.print(table)
[docs]@dataclass
class MomentCurvatureResults:
r"""Class for storing moment curvature results.
:param theta: Angle (in radians) the neutral axis makes with the horizontal
:param n_target: Target axial force
axis (:math:`-\pi \leq \theta \leq \pi`)
:param kappa: List of curvatures
:param n: List of axial forces
:param m_x: List of bending moments about the x-axis
:param m_y: List of bending moments about the y-axis
:param m_xy: List of resultant bending moments
:param failure_geometry: Geometry object of the region of the cross-section that
failed, ending the moment curvature analysis
:param convergence: The critical ratio between the strain and the failure strain
within the cross-section for each curvature step in the analysis. A value of one
indicates failure.
"""
# results
theta: float
n_target: float
kappa: List[float] = field(default_factory=list)
n: List[float] = field(default_factory=list)
m_x: List[float] = field(default_factory=list)
m_y: List[float] = field(default_factory=list)
m_xy: List[float] = field(default_factory=list)
failure_geometry: CPGeom = field(init=False, repr=False)
convergence: List[float] = field(default_factory=list)
# for analysis
_kappa: float = field(default=0, repr=False)
_n_i: float = field(default=0, repr=False)
_m_x_i: float = field(default=0, repr=False)
_m_y_i: float = field(default=0, repr=False)
_failure: bool = field(default=False, repr=False)
_failure_convergence: float = field(default=0, repr=False)
[docs] def plot_results(
self,
m_scale: float = 1e-6,
fmt: str = "o-",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots the moment curvature results.
:param m_scale: Scaling factor to apply to bending moment
:param fmt: Plot format string
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
# scale moments
moments = np.array(self.m_xy) * m_scale
# create plot and setup the plot
with plotting_context(title="Moment-Curvature", **kwargs) as (
fig,
ax,
):
ax.plot(self.kappa, moments, fmt) # type: ignore
plt.xlabel("Curvature")
plt.ylabel("Moment")
plt.grid(True)
return ax
[docs] @staticmethod
def plot_multiple_results(
moment_curvature_results: List[MomentCurvatureResults],
labels: List[str],
m_scale: float = 1e-6,
fmt: str = "o-",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots multiple moment curvature results.
:param moment_curvature_results: List of moment curvature results objects
:param labels: List of labels for each moment curvature diagram
:param m_scale: Scaling factor to apply to bending moment
:param fmt: Plot format string
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
# create plot and setup the plot
with plotting_context(title="Moment-Curvature", **kwargs) as (
fig,
ax,
):
idx = 0
# for each M-k curve
for idx, mk_result in enumerate(moment_curvature_results):
# scale results
kappas = np.array(mk_result.kappa)
moments = np.array(mk_result.m_xy) * m_scale
ax.plot(kappas, moments, fmt, label=labels[idx]) # type: ignore
plt.xlabel("Curvature")
plt.ylabel("Moment")
plt.grid(True)
# if there is more than one curve show legend
if idx > 0:
ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) # type: ignore
return ax
[docs] def plot_failure_geometry(
self,
title: str = "Failure Geometry",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots the geometry that fails in the moment curvature analysis.
:param title: Plot title
:param kwargs: Passed to
:meth:`~sectionproperties.pre.geometry.CompoundGeometry.plot_geometry`
:return: Matplotlib axes object
"""
return self.failure_geometry.plot_geometry(title=title, **kwargs)
[docs] def get_curvature(
self,
moment: float,
) -> float:
"""Given a moment, uses the moment-curvature results to interpolate a curvature.
:param moment: Bending moment at which to obtain curvature
:raises ValueError: If supplied moment is outside bounds of moment-curvature
results.
:return: Curvature
"""
# check moment is within bounds of results
m_min = min(self.m_xy)
m_max = max(self.m_xy)
if moment > m_max or moment < m_min:
raise ValueError(
"moment must be within the bounds of the moment-curvature results."
)
f_kappa = interp1d(
x=self.m_xy,
y=self.kappa,
kind="linear",
)
return float(f_kappa(moment))
[docs]@dataclass(order=True)
class UltimateBendingResults:
r"""Class for storing ultimate bending results.
:param theta: Angle (in radians) the neutral axis makes with the horizontal
axis (:math:`-\pi \leq \theta \leq \pi`)
:param d_n: Ultimate neutral axis depth
:param k_u: Neutral axis parameter *(d_n / d)*
:param n: Resultant axial force
:param m_x: Resultant bending moment about the x-axis
:param m_y: Resultant bending moment about the y-axis
:param m_xy: Resultant bending moment
:param label: Result label
"""
# bending angle
theta: float
# ultimate neutral axis depth
d_n: float = 0
k_u: float = 0
# resultant actions
n: float = 0
m_x: float = 0
m_y: float = 0
m_xy: float = 0
# label
label: Optional[str] = field(default=None, compare=False)
[docs] def print_results(
self,
fmt: str = "8.6e",
):
"""Prints the ultimate bending results to the terminal.
:param fmt: Number format
"""
table = Table(title="Ultimate Bending Results")
table.add_column("Property", justify="left", style="cyan", no_wrap=True)
table.add_column("Value", justify="right", style="green")
if self.label:
table.add_row("Label", self.label)
table.add_row("Bending Angle - theta", "{:>{fmt}}".format(self.theta, fmt=fmt))
table.add_row("Neutral Axis Depth - d_n", "{:>{fmt}}".format(self.d_n, fmt=fmt))
table.add_row(
"Neutral Axis Parameter - k_u", "{:>{fmt}}".format(self.k_u, fmt=fmt)
)
table.add_row("Axial Force", "{:>{fmt}}".format(self.n, fmt=fmt))
table.add_row("Bending Capacity - m_x", "{:>{fmt}}".format(self.m_x, fmt=fmt))
table.add_row("Bending Capacity - m_y", "{:>{fmt}}".format(self.m_y, fmt=fmt))
table.add_row("Bending Capacity - m_xy", "{:>{fmt}}".format(self.m_xy, fmt=fmt))
console = Console()
console.print(table)
[docs]@dataclass
class MomentInteractionResults:
"""Class for storing moment interaction results.
:param results: List of ultimate bending result objects
"""
results: List[UltimateBendingResults] = field(default_factory=list)
[docs] def sort_results(self) -> None:
"""Sorts the results by decreasing axial force."""
self.results.sort(reverse=True)
# remove duplicates from sorted list
new_results = []
for res in self.results:
if res not in new_results:
new_results.append(res)
self.results = new_results
[docs] def get_results_lists(
self,
moment: str,
) -> Tuple[List[float], List[float]]:
"""Returns a list of axial forces and moments.
:param moment: Which moment to plot, acceptable values are ``"m_x"``, ``"m_y"``
or ``"m_xy"``
:return: List of axial forces and moments *(n, m)*
"""
# build list of results
n_list = []
m_list = []
for result in self.results:
n_list.append(result.n)
if moment == "m_x":
m_list.append(result.m_x)
elif moment == "m_y":
m_list.append(result.m_y)
elif moment == "m_xy":
m_list.append(result.m_xy)
else:
raise ValueError(f"{moment} not an acceptable value for moment.")
return n_list, m_list
[docs] def plot_diagram(
self,
n_scale: float = 1e-3,
m_scale: float = 1e-6,
moment: str = "m_x",
fmt: str = "o-",
labels: bool = False,
label_offset: bool = False,
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots a moment interaction diagram.
:param n_scale: Scaling factor to apply to axial force
:param m_scale: Scaling factor to apply to the bending moment
:param moment: Which moment to plot, acceptable values are ``"m_x"``, ``"m_y"``
or ``"m_xy"``
:param fmt: Plot format string
:param labels: If set to True, also plots labels on the diagram
:param label_offset: If set to True, attempts to offset the label from the
diagram
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
# create plot and setup the plot
with plotting_context(title="Moment Interaction Diagram", **kwargs) as (
fig,
ax,
):
# get results
n_list, m_list = self.get_results_lists(moment=moment)
# scale results
forces = np.array(n_list) * n_scale
moments = np.array(m_list) * m_scale
# plot diagram
ax.plot(moments, forces, fmt) # type: ignore
# plot labels
if labels:
if label_offset:
# compute gradients of curve and aspect ratio of plot
grad = np.gradient([moments, forces], axis=1)
x_diff = ax.get_xlim() # type: ignore
y_diff = ax.get_ylim() # type: ignore
ar = (y_diff[1] - y_diff[0]) / (x_diff[1] - x_diff[0])
for idx, m in enumerate(m_list):
if self.results[idx].label:
# get x,y position on plot
x = m * m_scale
y = n_list[idx] * n_scale
if label_offset:
# calculate text offset
grad_pt = grad[1, idx] / grad[0, idx] / ar # type: ignore
if grad_pt == 0:
norm_angle = np.pi / 2
else:
norm_angle = np.arctan2(-1 / grad_pt, 1)
x_t = np.cos(norm_angle) * 20
y_t = np.sin(norm_angle) * 20
annotate_dict = {
"xytext": (x_t, y_t),
"textcoords": "offset points",
"arrowprops": dict(
arrowstyle="->",
connectionstyle="angle,angleA=0,angleB=90,rad=10",
),
"bbox": dict(boxstyle="round", fc="0.8"),
}
else:
annotate_dict = {}
# plot text
ax.annotate( # type: ignore
text=self.results[idx].label, xy=(x, y), **annotate_dict
)
plt.xlabel("Bending Moment")
plt.ylabel("Axial Force")
plt.grid(True)
return ax
[docs] @staticmethod
def plot_multiple_diagrams(
moment_interaction_results: List[MomentInteractionResults],
labels: List[str],
n_scale: float = 1e-3,
m_scale: float = 1e-6,
moment: str = "m_x",
fmt: str = "o-",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots multiple moment interaction diagrams.
:param moment_interaction_results: List of moment interaction results objects
:param labels: List of labels for each moment interaction diagram
:param n_scale: Scaling factor to apply to axial force
:param m_scale: Scaling factor to apply to bending moment
:param moment: Which moment to plot, acceptable values are ``"m_x"``, ``"m_y"``
or ``"m_xy"``
:param fmt: Plot format string
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
# create plot and setup the plot
with plotting_context(title="Moment Interaction Diagram", **kwargs) as (
fig,
ax,
):
idx = 0
# for each M-N curve
for idx, mi_result in enumerate(moment_interaction_results):
n_list, m_list = mi_result.get_results_lists(moment=moment)
# scale results
forces = np.array(n_list) * n_scale
moments = np.array(m_list) * m_scale
ax.plot(moments, forces, fmt, label=labels[idx]) # type: ignore
plt.xlabel("Bending Moment")
plt.ylabel("Axial Force")
plt.grid(True)
# if there is more than one curve show legend
if idx > 0:
ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) # type: ignore
return ax
[docs] def point_in_diagram(
self,
n: float,
m: float,
moment: str = "m_x",
) -> bool:
"""Determines whether or not the combination of axial force and moment lies
within the moment interaction diagram.
:param n: Axial force
:param m: Bending moment
:param moment: Which moment to analyse, acceptable values are ``"m_x"``,
``"m_y"`` or ``"m_xy"``
:returns: True, if combination of axial force and moment is within the diagram
"""
# get results
n_list, m_list = self.get_results_lists(moment=moment)
# create a polygon from points on diagram
poly_points = []
for idx, mom in enumerate(m_list):
poly_points.append((mom, n_list[idx]))
poly = Polygon(poly_points)
point = Point(m, n)
return poly.contains(point)
[docs]@dataclass
class BiaxialBendingResults:
"""Class for storing biaxial bending results.
:param n: Net axial force
:param results: List of ultimate bending result objects
"""
n: float
results: List[UltimateBendingResults] = field(default_factory=list)
[docs] def get_results_lists(
self,
) -> Tuple[List[float], List[float]]:
"""Returns a list and moments about the ``x`` and ``y`` axes.
:return: List of axial forces and moments *(mx, my)*
"""
# build list of results
m_x_list = []
m_y_list = []
for result in self.results:
m_x_list.append(result.m_x)
m_y_list.append(result.m_y)
return m_x_list, m_y_list
[docs] def plot_diagram(
self,
m_scale: float = 1e-6,
fmt: str = "o-",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots a biaxial bending diagram.
:param m_scale: Scaling factor to apply to bending moment
:param fmt: Plot format string
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
m_x_list, m_y_list = self.get_results_lists()
# create plot and setup the plot
with plotting_context(
title=f"Biaxial Bending Diagram, $N = {self.n:.3e}$", **kwargs
) as (
fig,
ax,
):
# scale results
m_x = np.array(m_x_list) * m_scale
m_y = np.array(m_y_list) * m_scale
ax.plot(m_x, m_y, fmt) # type: ignore
plt.xlabel("Bending Moment $M_x$")
plt.ylabel("Bending Moment $M_y$")
plt.grid(True)
return ax
[docs] @staticmethod
def plot_multiple_diagrams_2d(
biaxial_bending_results: List[BiaxialBendingResults],
labels: Optional[List[str]] = None,
m_scale: float = 1e-6,
fmt: str = "o-",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots multiple biaxial bending diagrams in a 2D plot.
:param biaxial_bending_results: List of biaxial bending results objects
:param labels: List of labels for each biaxial bending diagram, if not provided
labels are axial forces
:param m_scale: Scaling factor to apply to bending moment
:param fmt: Plot format string
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
# create plot and setup the plot
with plotting_context(title="Biaxial Bending Diagram", **kwargs) as (
fig,
ax,
):
idx = 0
# generate default labels
if labels is None:
labels = []
default_labels = True
else:
default_labels = False
# for each M-N curve
for idx, bb_result in enumerate(biaxial_bending_results):
m_x_list, m_y_list = bb_result.get_results_lists()
# scale results
m_x_list = np.array(m_x_list) * m_scale
m_y_list = np.array(m_y_list) * m_scale
# generate default labels
if default_labels:
labels.append(f"N = {bb_result.n:.3e}")
ax.plot(m_x_list, m_y_list, fmt, label=labels[idx]) # type: ignore
plt.xlabel("Bending Moment $M_x$")
plt.ylabel("Bending Moment $M_y$")
plt.grid(True)
# if there is more than one curve show legend
if idx > 0:
ax.legend(loc="center left", bbox_to_anchor=(1, 0.5)) # type: ignore
return ax
[docs] @staticmethod
def plot_multiple_diagrams_3d(
biaxial_bending_results: List[BiaxialBendingResults],
n_scale: float = 1e-3,
m_scale: float = 1e-6,
fmt: str = "-",
) -> matplotlib.axes.Axes: # type: ignore
"""Plots multiple biaxial bending diagrams in a 3D plot.
:param biaxial_bending_results: List of biaxial bending results objects
:param n_scale: Scaling factor to apply to axial force
:param m_scale: Scaling factor to apply to bending moment
:param fmt: Plot format string
:return: Matplotlib axes object
"""
# make 3d plot
fig = plt.figure()
ax = plt.axes(projection="3d")
# for each curve
for bb_result in biaxial_bending_results:
m_x_list, m_y_list = bb_result.get_results_lists()
# scale results
n_list = bb_result.n * n_scale * np.ones(len(m_x_list))
m_x_list = np.array(m_x_list) * m_scale
m_y_list = np.array(m_y_list) * m_scale
ax.plot3D(m_x_list, m_y_list, n_list, fmt) # type: ignore
ax.set_xlabel("Bending Moment $M_x$")
ax.set_ylabel("Bending Moment $M_y$")
ax.set_zlabel("Axial Force $N$") # type: ignore
plt.show()
return ax
[docs] def point_in_diagram(
self,
m_x: float,
m_y: float,
) -> bool:
"""Determines whether or not the combination of bending moments lies within the
biaxial bending diagram.
:param m_x: Bending moment about the x-axis
:param m_y: Bending moment about the y-axis
:returns: True, if combination of bendings moments is within the diagram
"""
# create a polygon from points on diagram
poly_points = []
for ult_res in self.results:
poly_points.append((ult_res.m_x, ult_res.m_y))
poly = Polygon(poly_points)
point = Point(m_x, m_y)
return poly.contains(point)
[docs]@dataclass
class StressResult:
"""Class for storing stress results.
The lever arm is computed to the elastic centroid.
:param concrete_analysis_sections: List of concrete analysis section objects
present in the stress analysis, which can be visualised by calling the
:meth:`~concreteproperties.analysis_section.AnalysisSection.plot_mesh` or
:meth:`~concreteproperties.analysis_section.AnalysisSection.plot_shape`
:param concrete_stresses: List of concrete stresses at the nodes of each concrete
analysis section
:param concrete_forces: List of net forces for each concrete analysis section and
its lever arm (``force``, ``d_x``, ``d_y``)
:param meshed_reinforcement_sections: List of meshed reinforcement section objects
present in the stress analysis
:param meshed_reinforcement_stresses: List of meshed reinforcement stresses at the
nodes of each meshed reinforcement analysis section
:param meshed_reinforcement_forces: List of net forces for each meshed reinforcement
analysis section and its lever arm (``force``, ``d_x``, ``d_y``)
:param lumped_reinforcement_geometries: List of lumped reinforcement geometry
objects present in the stress analysis
:param lumped_reinforcement_stresses: List of lumped reinforcement stresses for
each lumped geometry
:param lumped_reinforcement_strains: List of lumped reinforcement strains for each
lumped geometry
:param lumped_reinforcement_forces: List of net forces for each lumped reinforcement
geometry and its lever arm (``force``, ``d_x``, ``d_y``)
:param strand_geometries: List of strand geometry objects present in the stress
analysis
:param strand_stresses: List of strand stresses for each strand
:param strand_strains: List of strand strains for each strand
:param strand_forces: List of net forces for each strand geometry and its lever arm
(``force``, ``d_x``, ``d_y``)
"""
concrete_section: ConcreteSection
concrete_analysis_sections: List[AnalysisSection]
concrete_stresses: List[np.ndarray]
concrete_forces: List[Tuple[float, float, float]]
meshed_reinforcement_sections: List[AnalysisSection]
meshed_reinforcement_stresses: List[np.ndarray]
meshed_reinforcement_forces: List[Tuple[float, float, float]]
lumped_reinforcement_geometries: List[CPGeom]
lumped_reinforcement_stresses: List[float]
lumped_reinforcement_strains: List[float]
lumped_reinforcement_forces: List[Tuple[float, float, float]]
strand_geometries: List[CPGeom] = field(default_factory=list)
strand_stresses: List[float] = field(default_factory=list)
strand_strains: List[float] = field(default_factory=list)
strand_forces: List[Tuple[float, float, float]] = field(default_factory=list)
_m_net: Optional[float] = field(default=None, repr=False)
[docs] def plot_stress(
self,
title: str = "Stress",
conc_cmap: str = "RdGy",
reinf_cmap: str = "bwr",
**kwargs,
) -> matplotlib.axes.Axes: # type: ignore
"""Plots concrete and steel stresses on a concrete section.
:param title: Plot title
:param conc_cmap: Colour map for the concrete stress
:param reinf_cmap: Colour map for the reinforcement stress
:param kwargs: Passed to :func:`~concreteproperties.post.plotting_context`
:return: Matplotlib axes object
"""
with plotting_context(
title=title,
aspect=True,
**dict(
kwargs, nrows=1, ncols=3, gridspec_kw={"width_ratios": [1, 0.08, 0.08]}
),
) as (fig, ax):
# plot background
self.concrete_section.plot_section(
background=True, **dict(kwargs, ax=fig.axes[0])
)
# set up the colormaps
cmap_conc = cm.get_cmap(name=conc_cmap)
cmap_reinf = cm.get_cmap(name=reinf_cmap)
# determine minimum and maximum stress values for the contour list
# add tolerance for plotting stress blocks
conc_sig_min = min([min(x) for x in self.concrete_stresses]) - 1e-12
conc_sig_max = max([max(x) for x in self.concrete_stresses]) + 1e-12
# if there is meshed reinforcement, calculate min and max
if self.meshed_reinforcement_stresses:
meshed_reinf_sig_min = (
min([min(x) for x in self.meshed_reinforcement_stresses]) - 1e-12
)
meshed_reinf_sig_max = (
max([max(x) for x in self.meshed_reinforcement_stresses]) + 1e-12
)
else:
meshed_reinf_sig_min = None
meshed_reinf_sig_max = None
# if there is lumped reinforcement, calculate min and max
if self.lumped_reinforcement_stresses or self.strand_stresses:
lumped_reinf_sig_min = min(
self.lumped_reinforcement_stresses + self.strand_stresses
)
lumped_reinf_sig_max = max(
self.lumped_reinforcement_stresses + self.strand_stresses
)
else:
lumped_reinf_sig_min = None
lumped_reinf_sig_max = None
# determine min and max reinforcement stresess
if (
meshed_reinf_sig_min
and meshed_reinf_sig_max
and lumped_reinf_sig_min
and lumped_reinf_sig_max
):
reinf_sig_min = min(meshed_reinf_sig_min, lumped_reinf_sig_min)
reinf_sig_max = max(meshed_reinf_sig_max, lumped_reinf_sig_max)
elif meshed_reinf_sig_min and meshed_reinf_sig_max:
reinf_sig_min = meshed_reinf_sig_min
reinf_sig_max = meshed_reinf_sig_max
elif lumped_reinf_sig_min and lumped_reinf_sig_max:
reinf_sig_min = lumped_reinf_sig_min
reinf_sig_max = lumped_reinf_sig_max
else:
reinf_sig_min = 0
reinf_sig_max = 0
# set up ticks
v_conc = np.linspace(conc_sig_min, conc_sig_max, 15, endpoint=True)
v_reinf = np.linspace(reinf_sig_min, reinf_sig_max, 15, endpoint=True)
if np.isclose(v_conc[0], v_conc[-1], atol=1e-12):
v_conc = 15
ticks_conc = None
else:
ticks_conc = v_conc
if np.isclose(v_reinf[0], v_reinf[-1], atol=1e-12):
ticks_reinf = None
reinf_tick_same = True
else:
ticks_reinf = v_reinf
reinf_tick_same = False
# plot the concrete stresses
for idx, sig in enumerate(self.concrete_stresses):
# check region has a force
if abs(self.concrete_forces[idx][0]) > 1e-8:
# create triangulation
triang_conc = tri.Triangulation(
self.concrete_analysis_sections[idx].mesh_nodes[:, 0],
self.concrete_analysis_sections[idx].mesh_nodes[:, 1],
self.concrete_analysis_sections[idx].mesh_elements[:, 0:3], # type: ignore
)
# plot the filled contour
trictr_conc = fig.axes[0].tricontourf(
triang_conc, sig, v_conc, cmap=cmap_conc, norm=CenteredNorm()
) # type: ignore
# plot a zero stress contour, supressing warning
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message="No contour levels were found within the data range.",
)
# set zero stress for neutral axis contour
zero_level = 0
if min(sig) > 0:
if min(sig) < 1e-3:
zero_level = min(sig) + 1e-12
if max(sig) < 0:
if max(sig) > -1e-3:
zero_level = max(sig) - 1e-12
if min(sig) == 0:
zero_level = 1e-12
if max(sig) == 0:
zero_level = -1e-12
CS = fig.axes[0].tricontour(
triang_conc,
sig,
[zero_level],
linewidths=1,
linestyles="dashed",
) # type: ignore
# plot the meshed reinforcement stresses
trictr_reinf = None
for idx, sig in enumerate(self.meshed_reinforcement_stresses):
# check region has a force
if abs(self.meshed_reinforcement_forces[idx][0]) > 1e-8:
# create triangulation
triang_reinf = tri.Triangulation(
self.meshed_reinforcement_sections[idx].mesh_nodes[:, 0],
self.meshed_reinforcement_sections[idx].mesh_nodes[:, 1],
self.meshed_reinforcement_sections[idx].mesh_elements[:, 0:3], # type: ignore
)
# plot the filled contour
trictr_reinf = fig.axes[0].tricontourf(
triang_reinf, sig, v_reinf, cmap=cmap_reinf, norm=CenteredNorm()
) # type: ignore
# plot a zero stress contour, supressing warning
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message="No contour levels were found within the data range.",
)
# set zero stress for neutral axis contour
zero_level = 0
if min(sig) > 0:
if min(sig) < 1e-3:
zero_level = min(sig) + 1e-12
if max(sig) < 0:
if max(sig) > -1e-3:
zero_level = max(sig) - 1e-12
if min(sig) == 0:
zero_level = 1e-12
if max(sig) == 0:
zero_level = -1e-12
CS = fig.axes[0].tricontour(
triang_reinf,
sig,
[zero_level],
linewidths=1,
linestyles="dashed",
) # type: ignore
# plot the lumped reinforcement stresses
lumped_reinf_patches = []
colours = []
for idx, sig in enumerate(self.lumped_reinforcement_stresses):
lumped_reinf_patches.append(
mpatches.Polygon(
xy=list(self.lumped_reinforcement_geometries[idx].geom.exterior.coords) # type: ignore
)
)
colours.append(sig)
for idx, sig in enumerate(self.strand_stresses):
lumped_reinf_patches.append(
mpatches.Polygon(
xy=list(self.strand_geometries[idx].geom.exterior.coords) # type: ignore
)
)
colours.append(sig)
patch = PatchCollection(lumped_reinf_patches, cmap=cmap_reinf)
patch.set_array(colours)
if reinf_tick_same:
patch.set_clim(vmin=0.99 * v_reinf[0], vmax=1.01 * v_reinf[-1])
else:
patch.set_clim(vmin=v_reinf[0], vmax=v_reinf[-1])
fig.axes[0].add_collection(patch) # type: ignore
# add the colour bars
fig.colorbar(
trictr_conc, # type: ignore
label="Concrete Stress",
format="%.2e",
ticks=ticks_conc,
cax=fig.axes[1],
)
if trictr_reinf:
mappable = trictr_reinf
else:
mappable = patch
fig.colorbar(
mappable,
label="Reinforcement Stress",
format="%.2e",
ticks=ticks_reinf,
cax=fig.axes[2],
)
return ax
[docs] def sum_forces(self) -> float:
"""Returns the sum of the internal forces.
:return: Sum of internal forces
"""
force_sum = 0
# sum concrete forces
for conc_force in self.concrete_forces:
force_sum += conc_force[0]
# sum meshed reinf forces
for meshed_reinf_force in self.meshed_reinforcement_forces:
force_sum += meshed_reinf_force[0]
# sum lumped reinf forces
for lumped_reinf_force in self.lumped_reinforcement_forces:
force_sum += lumped_reinf_force[0]
# sum strand forces
for strand_force in self.strand_forces:
force_sum += strand_force[0]
return force_sum
[docs] def sum_moments(self) -> Tuple[float, float, float]:
"""Returns the sum of the internal moments.
:return: Sum of internal moments about each axis and resultant moment
(``m_x``, ``m_y``, ``m``)
"""
moment_sum_x = 0
moment_sum_y = 0
# sum concrete moments
for conc_force in self.concrete_forces:
moment_sum_x += conc_force[0] * conc_force[2]
moment_sum_y += conc_force[0] * conc_force[1]
# sum meshed reinf moments
for meshed_reinf_force in self.meshed_reinforcement_forces:
moment_sum_x += meshed_reinf_force[0] * meshed_reinf_force[2]
moment_sum_y += meshed_reinf_force[0] * meshed_reinf_force[1]
# sum lumped reinf moments
for lumped_reinf_force in self.lumped_reinforcement_forces:
moment_sum_x += lumped_reinf_force[0] * lumped_reinf_force[2]
moment_sum_y += lumped_reinf_force[0] * lumped_reinf_force[1]
# sum strand moments
for strand_force in self.strand_forces:
moment_sum_x += strand_force[0] * strand_force[2]
moment_sum_y += strand_force[0] * strand_force[1]
moment_sum = np.sqrt(moment_sum_x * moment_sum_x + moment_sum_y * moment_sum_y)
return moment_sum_x, moment_sum_y, moment_sum
[docs] def get_concrete_stress_limits(self) -> Tuple[float, float]:
"""Returns the minimum and maximum concrete stress.
:return: Minimum concrete stress, maximum concrete stress
"""
min_stress = 0
max_stress = 0
for idx, stress_list in enumerate(self.concrete_stresses):
if idx == 0:
min_stress = stress_list.min()
max_stress = stress_list.max()
else:
min_stress = min(min_stress, stress_list.min())
max_stress = max(max_stress, stress_list.max())
return min_stress, max_stress