diff --git a/.github/workflows/coveralls.yml b/.github/workflows/coveralls.yml index d76e4e6..c6d7cd6 100644 --- a/.github/workflows/coveralls.yml +++ b/.github/workflows/coveralls.yml @@ -29,7 +29,7 @@ jobs: - name: Test with pytest run: | - python -m pytest --cov=efficalc --cov-report=xml tests + python -m pytest tests --cov=efficalc --cov-report=xml tests - name: Coveralls uses: coverallsapp/github-action@v2 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index dc48cd3..eb34164 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -37,4 +37,4 @@ jobs: flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - name: Test with pytest run: | - python -m pytest + python -m pytest tests diff --git a/examples/conc_col_pmm/README.md b/examples/conc_col_pmm/README.md new file mode 100644 index 0000000..74c3640 --- /dev/null +++ b/examples/conc_col_pmm/README.md @@ -0,0 +1,281 @@ +# Scope +This program supports the analysis of rectangular, symmetric, perimeter-reinforced concrete columns in accordance with ACI 318-19 and using the exact capacity method. + +# Run Guide + +## Getting Started +The `conc_col_pmm` folder is meant to be a standalone design tool. To get set up for running it on your local machine, follow these steps: + +1. Download or copy the entire `conc_col_pmm` folder into your desired working directory. +2. Open a command line shell in the working directory. Make sure you have python installed and accessible. +3. Initialize a virtual environment (e.g. `python -m venv .venv && .venv/Scripts/activate.ps1 ` ) +4. Install requirements (e.g. `pip install -r requirements.txt`) +5. Run tests to make sure setup is complete: `pytest conc_col_pmm/tests` + +## Designing a column +The `conc_col_pmm/tests/visual_tests` folder has examples for running various parts of the concrete column tool including + +* `visual_test_document_wrapper.py` for viewing a complete calculation report +* `visual_test_pmm_plotter_plotly.py` for viewing a 3D PMM plot +* `visual_test_point_plotter.py` for viewing 2D PM plots for specific load cases + +These example files can be run and viewed with `python -m conc_col_pmm.tests.visual_tests.`. For example: +`python -m conc_col_pmm.tests.visual_tests.visual_test_document_wrapper` + + +# Program Features +- **PMM Diagrams**: Creates 3D Axial-Moment-Moment interaction diagrams. +- **PM Diagrams**: For each load case entered, creates a 2D PM interaction diagram including the location of the load case on the PM axes relative to the capacity curve. +- **Demand-to-Capacity Ratios (DCRs)**: Calculates the PM Vector DCR (defined below) for each load case entered. +- **Calculation Reports**: Generates a calc report for the column showing the full capacity calculation for all selected load cases, including cross-section diagrams, code references, and DCR calculations. + +# Definitions and Theory +- **Axial force ($P$)**: The force in the column parallel to its length, where a positive value indicates compression. +- **Ultimate axial force ($P_u$)**: The axial force applied to a column (typically calculated using a factored load combination). +- **Nominal axial capacity ($P_n$)**: The calculated axial load capacity, without a safety factor. +- **Moments ($M_x$, $M_y$)**: These measure moments about axes indicated by their subscripts, i.e., positive $M_{x}$ indicates compression in the +y region and positive $M_{y}$ indicates compression in the +x region. +- **Ultimate moments ($M_{ux}$, $M_{uy}$)**: The moments applied to a column (typically calculated using a factored load combination). +- **Nominal moment capacity ($M_{nx}$, $M_{ny}$)**: The calculated moment capacities, without a safety factor. +- **Load case**: A given combination of axial load and moments ($P_{u}$, $M_{ux}$, $M_{uy}$). +- **Eccentricity ($e_x$, $e_y$)**: This is the location that a load would have to have on the x or y axis to produce an equivalent moment to the moment at the given load case with the same axial load. + +$$ +e_x=\frac{M_x}{P} +$$ + +$$ +e_y=\frac{M_y}{P} +$$ +- **Eccentricity angle (λ)**: The angle taken clockwise from the positive y-axis to a vector from the column centroid to the location where an equivalent axial load could be applied that would be equivalent to the combination of $P$, $M_{x}$, and $M_{y}$ for a given load case. + +$$ +\lambda = \arctan\left(\frac{e_x}{e_y}\right) = \arctan\left(\frac{M_{ny}}{M_{nx}}\right) +$$ +- **Neutral axis**: The line across the column section on which strain is assumed to be 0. +- **Neutral axis angle (θ)**: The angle from the positive x-axis to the neutral axis. Counter-clockwise is taken as positive, and the region of compression is on the side of the neutral axis further counter-clockwise than θ. The neutral axis angle can be understood in the graphic below, which shows the neutral axis as a dashed line and the equivalent compression zone as gray. + + +- **Neutral axis depth (c)**: The distance from the neutral axis to a parallel line passing through the column corner in maximum compression. +- **PMM diagram**: The 3D plot of the PMM surface, which describes the capacity of a given column in $P$, $M_x$, $M_y$ (any load case whose plot falls inside the PMM surface is within capacity, and any load case falling outside the PMM surface exceeds capacity). +- **PM diagram**: The 2D plot of a cut of the PMM surface at a given eccentricity angle (λ). This diagram is useful for visualizing the capacity of the column relative to demand for a given load case. +- **PM Vector DCR**: The ratio of the length of the demand vector (in PMM space) to the length of a parallel vector beginning at the origin and continuing until it reaches the capacity surface. +- **Axial to Moment Angle (α)**: This is a custom-defined variable used in this program that describes the height of a load case above the $M_x-M_y$ plane. It is a proxy for the inverse of the eccentricity resultant. + +$$ +\alpha = \arctan\left(\frac{P_n}{\sqrt{M_{nx}^2+M_{ny}^2}}\right) +$$ + +It seems intuitive that the neutral axis should be parallel to the axis of the resultant moment, which would mean the relation λ=-θ would hold. However, this holds only in special cases, which means determining the neutral axis angle required to produce a given eccentricity is not straightforward. For more information see [1]. + +# How It Works +- **PMM Diagrams**: The PMM diagrams created by this program compose a mesh of capacity points which are evenly spaced in both the vertical load ($P$) direction and in their angles about the origin (λ). To achieve this even spacing, it is necessary to find points on the PMM surface that have a given combination of λ and $P$. $P$ tends to increase as the neutral axis depth increases and λ tends to increase as θ decreases, but neither of the output variables (λ and $P$) can be calculated in closed form. This means that the domain of the two input variables (θ and c) must be searched to find the target point on the PMM diagram. The search algorithm is [below](#search-algorithm). +- **PM Diagrams**: This program uses sets of control points interpolated from the points on the PMM diagram to create PM diagrams for each load case. +- **DCRs**: These are calculated by finding a point on the PMM surface such that a vector from the origin to that point is parallel to a vector from the origin to the PMM point for the given load case and then taking a ratio of the lengths of the two vectors. To find the capacity point, a search of the PMM surface is again required, but in this case, the target variables are λ and α. Searching this domain is equivalent to searching in the domain of spherical coordinates. Unlike in the case of the point search for the PMM diagram, this search is performed on the fully-factored PMM surface (including the plateau). The search algorithm is [below](#search-algorithm). +- **Calculation Reports**: The program uses the efficalc library to generate calc reports. + +# Search Algorithm +## Summary +In both search problems, there are two input and two output variables and the target output variables are known. The algorithm is given a starting point, and it proceeds as follows: +1. Calculate first derivatives—the full 4x4 Jacobian—at the current input point using finite differences. +2. Use the derivatives as linear approximations for the two output variables as functions of the two variables and solve for the input point at which both output variables are expected to equal their target values. +3. Move to the input point calculated in (2) and repeat from (1). + +## Update Method +For simplicity, assume that the two input variables are x and y and the two output functions are f and g, where f and g have been shifted so that the target outputs are f=0 and g=0. Then the Jacobian is as follows: + +$$ +J = +\begin{bmatrix} +\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ +\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} +\end{bmatrix} +$$ + +Then the linear approximator can be written: + +$$ +\begin{bmatrix} +f_1 \\ +g_1 +\end{bmatrix}= +\begin{bmatrix} +\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ +\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} +\end{bmatrix} +\begin{bmatrix} +x \\ +y +\end{bmatrix}+ +\begin{bmatrix} +f_0 \\ +g_0 +\end{bmatrix} +$$ + +Since the target is f=0, g=0, this is equivalent to the following system: + +$$ +-\begin{bmatrix} +f_0 \\ +g_0 +\end{bmatrix}= +\begin{bmatrix} +\frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ +\frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} +\end{bmatrix} +\begin{bmatrix} +x \\ +y +\end{bmatrix} +$$ + +This equation yields solutions for x and y: + +$$ +x = -\frac{\frac{\partial g}{\partial y} f_0 - \frac{\partial f}{\partial y} g_0}{\frac{\partial f}{\partial x} \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \frac{\partial g}{\partial x}} +$$ + +$$ +y = -\frac{-\frac{\partial g}{\partial x} f_0 + \frac{\partial f}{\partial x} g_0}{\frac{\partial f}{\partial x} \frac{\partial g}{\partial y} - \frac{\partial f}{\partial y} \frac{\partial g}{\partial x}} +$$ + +This program uses the formulas above to calculate the next guess of both inputs at each iteration. The situation can be visualized with the following plot of the projected zero-contours of the two functions f and g, where both zero-contours are estimated by calculating the gradient at the current point. The solution to the linear system above effectively finds the intersection between the two zero-contours, which is the target point. + + + + +# Summary of Program Structure by Package and Module + +## 1. `calc_document` +Contains modules linked to the `efficalc` package which are used for generating the calc report. + +### 1.1. `add_col_inputs_document` +Adds the column inputs and assumptions to the calc report. + +### 1.2. `col_inputs` +Collects information from the user about the column and loads. Calls `full_calc_document` to begin creating the calc report. + +### 1.3. `dcr_calc_runner` +Calculates DCRs for all load cases and adds calculations to the calc report for applicable load cases. + +### 1.4. `document_wrapper` +Creates the calc report. It calls `col_inputs`, and from there, all information is added to the calc report. + +### 1.5. `full_calc_document` +Receives a column and load cases, then runs the calculations for the column capacity and DCRs for all load cases. + +### 1.6. `results_summary` +Creates a table showing the DCRs for all load cases and shows the max DCR. + +### 1.7. `show_dcr_calc` +Adds the calculation of a particular DCR to the calc report. Optionally called depending on whether the user selects a given load case to be shown. + +### 1.8. `try_axis_document` +Adds the calculations for the reaction of a column to bending on a given neutral axis to the calc report. + +### 1.9. `plotting` (sub-package) +Contains plotting functions. + +#### 1.9.1. `get_capacity` +Accepts parameters like the quarter PMM mesh and a loading point, then returns a list of resultant moment and axial points which form the PM diagram at the angle of the given load point. + +#### 1.9.2. `get_pmm_data` +Creates an instance of the PMM class containing the data for a given column's PMM diagram. + +#### 1.9.3. `PMM` +Defines a dataclass for storing the data needed for plotting the PMM diagram. + +#### 1.9.4. `pmm_mesh` +Creates the mesh for the PMM diagram by iterating over the range of axial loads and λ. + +#### 1.9.5. `pmm_plotter_plotly` +Creates a Plotly figure for the column’s PMM diagram. + +#### 1.9.6. `point_plotter` +Creates a Matplotlib figure showing the PM diagram for a given load case, along with the point for the given load case. + +#### 1.9.7. `pure_mx_my_plotter` +Adds PM diagrams to the calc report showing cuts of the PMM diagram that align with both the x and y axes (indicating moment purely about the x and y axes). + +## 2. `col` +Contains functions related to the definition of a given concrete column. + +### 2.1. `assign_max_min` +Performs calculations for the maximum and minimum axial capacity of a given column, adds them to the calc report, and assigns the calculated values to the given `Column` object. + +### 2.2. `column` +Contains the class defining a `Column` object, including various properties. + +### 2.3. `col_canvas` (sub-package) +Contains plotting functions. + +#### 2.3.1. `draw_column_comp_zone` +Draws the cross-section of the column with the full equivalent compression zone labeled. + +#### 2.3.2. `draw_column_with dimensions` +Draws the cross-section of the column with dimensions and rebar information shown. + +#### 2.3.3. `draw_column_with_triangle` +Draws the cross-section of the column with a triangular compression area labeled. + +#### 2.3.4. `draw_plain_column` +Draws the cross-section of the column on an `efficalc` Canvas, including rebar. + +## 3. `constants` +Contains constants used by other packages. + +### 3.1. `rebar_data` +Stores data on rebar properties. + +## 4. `pmm_search` +Contains functions used for searching the PMM surface for target points. + +### 4.1. `ecc_search` +Used to find a PMM point for DCR calculation. + +#### 4.1.1. `change_ecc` +Calculates the next iteration point in the search. + +#### 4.1.2. `get_dcr_ecc` +Runs the point search for finding the capacity point for a given load point and calculates the DCR. Optionally adds the DCR calculation to the calc report. + +#### 4.1.3. `get_error_ecc` +Calculates the distance of the current iteration point from the target point in normalized λ-α space. + +#### 4.1.4. `limit_comp_ecc` +Calculates the results for a given iteration point by running `try_axis`. Also limits the axial load to a range in which nonzero derivatives can be calculated. + +#### 4.1.5. `point_search_ecc` +Controls the convergence of the search for a point for a DCR calculation. + +### 4.2. `load_search` +Used to find a PMM point for the PMM diagram. + +#### 4.2.1. `bisect_load` +Searches for a point aligned with the Mx or My axes by changing only the neutral axis depth while holding λ constant. + +#### 4.2.2. `change_load` +Calculates the next iteration point in the search. + +#### 4.2.3. `get_error_load` +Calculates the distance of the current iteration point from the target point in normalized λ-P space. + +#### 4.2.4. `limit_comp_load` +Calculates the results for a given iteration point by running `try_axis`. Also limits the axial load to a range in which nonzero derivatives can be calculated. + +#### 4.2.5. `point_search_load` +Controls the convergence of the search for a point for creating a PMM diagram. + +#### 4.2.6. `starting_pts` +Selects starting points for the `bisect_load` algorithm. Both starting points are near the initial guess point. + +## 5. `struct_analysis` + +### 5.1. `triangles` +Contains functions for calculating the area and centroid of a triangle, for use in `try_axis`. + +### 5.2. `try_axis` +Calculates the reaction of a column (axial load and moments) due to bending on a given neutral axis angle and depth. + +# Reference +[1] Design of Concrete Structures, 15th ed. Darwin, Dolan, and Nilson. McGraw, 2016. diff --git a/examples/conc_col_pmm/__init__.py b/examples/conc_col_pmm/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/calc_document/__init__.py b/examples/conc_col_pmm/calc_document/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/calc_document/calculation.py b/examples/conc_col_pmm/calc_document/calculation.py new file mode 100644 index 0000000..e8f9dfd --- /dev/null +++ b/examples/conc_col_pmm/calc_document/calculation.py @@ -0,0 +1,157 @@ +from efficalc import Assumption, Calculation, Heading, Input, InputTable, Title + +from ..col.column import Column +from ..constants.concrete_data import MAX_CONCRETE_STRAIN +from ..constants.rebar_data import REBAR_SIZES, REBAR_STRENGTHS, STEEL_E, rebar_area +from ..pmm_search.load_combo import LoadCombination, is_yes +from .column_inputs import ColumnInputs +from .full_calc_document import calculation as full_calc + +# TODO: this should return all info needed to plot visual tests, possibly take input params for defaults + + +# this function accepts inputs from the user and passes them to "full_calc_document" +def calculation( + default_loads: list[list] = [[3000, -200, 100, "yes"]], col=ColumnInputs() +): + Title("Concrete Column Biaxial Bending Calculation Report") + + Heading("Column Inputs") + w = Input("w", col.w, "in", description="Column section width (x dimension)") + h = Input("h", col.h, "in", description="Column section height (y dimension)") + + # zero spaces + bar_size = Input( + "", + col.bar_size, + "", + description="Longitudinal rebar size (Imperial)", + input_type="select", + select_options=REBAR_SIZES, + ) + + # one space + bar_cover = Input(" ", col.bar_cover, "in", description="Longitudinal rebar cover") + + # TODO: should clear cover account for the shear reinforcement? + cover_options = [ + "Center", + "Edge", + ] + + # 2 spaces + cover_type = Input( + " ", + default_value=cover_options[0] if col.cover_to_center else cover_options[1], + unit="", + description="Cover is to bar center or bar edge (clear cover)", + input_type="select", + select_options=cover_options, + ) + transverse_options = ["Spiral", "Tied"] + # 3 spaces + transverse_type = Input( + " ", + default_value=( + transverse_options[0] if col.spiral_reinf else transverse_options[1] + ), + unit="", + description="Transverse reinforcement type", + input_type="select", + select_options=transverse_options, + ) + + # 4 spaces + bars_x = Input( + "n_x", + col.bars_x, + "", + description="Number of bars on the top/bottom edges", + num_step=1, + ) + + # 5 spaces + bars_y = Input( + "n_y", + col.bars_y, + "", + description="Number of bars on the left/right edges", + num_step=1, + ) + fc = Input("f^{\prime}_c", 8000, "psi", description="Concrete strength") + + fy = Input( + "f_y", + col.fy, + "ksi", + description="Steel strength", + input_type="select", + select_options=REBAR_STRENGTHS, + ) + + headers = [ + "Pu (kip)", + "Mux (kip-ft)", + "Muy (kip-ft)", + "Show in Calc Report (yes/no)", + ] + load_table = InputTable( + default_loads, headers, "Load Cases", False, False, numbered_rows=True + ) + + load_combos = [ + LoadCombination(idx + 1, load[0], load[1], load[2], is_yes(load[3])) + for idx, load in enumerate(load_table.data) + ] + + # above were the efficalc Inputs from the user, and below, some additional inputs and assumptions + # are created and added to the calc report + + A_b = Calculation( + "A_{\\mathrm{bar}}", + rebar_area(bar_size.get_value()), + "in^2", + description="Area of one bar", + ) + + E_s = Calculation( + "E_s", STEEL_E, "ksi", "Steel modulus of elasticity", "ACI 318-19 20.2.2.2" + ) + e_c = Calculation( + "\\epsilon_u", + MAX_CONCRETE_STRAIN, + "", + "Concrete strain at f'c", + "ACI 318-19 22.2.2.1", + ) + + Heading("Assumptions") + Assumption("ACI 318-19 controls the design") + Assumption("Reinforcement is non-prestressed") + Assumption( + "Lap splices of longitudinal reinforcement are in accordance with ACI 318-19 Table 10.7.5.2.2" + ) + Assumption( + "Strain in concrete and reinforcement is proportional to distance from the neutral axis, per ACI 318-19 22.2.1.2 " + ) + + cover_to_center = cover_type == "Center" + spiral_reinf = transverse_type == "Spiral" + + column = Column( + w, + h, + bar_size, + bar_cover, + bars_x, + bars_y, + fc, + fy, + cover_to_center, + spiral_reinf, + A_b, + E_s, + e_c, + ) + + full_calc(column, load_combos) diff --git a/examples/conc_col_pmm/calc_document/column_capacities.py b/examples/conc_col_pmm/calc_document/column_capacities.py new file mode 100644 index 0000000..f4d5084 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/column_capacities.py @@ -0,0 +1,12 @@ +import dataclasses + +from latexexpr_efficalc import Variable + +from efficalc import Calculation + + +@dataclasses.dataclass +class ColumnCapacities: + Mx: Calculation | Variable + My: Calculation | Variable + P: Calculation | Variable diff --git a/examples/conc_col_pmm/calc_document/column_inputs.py b/examples/conc_col_pmm/calc_document/column_inputs.py new file mode 100644 index 0000000..3c76faa --- /dev/null +++ b/examples/conc_col_pmm/calc_document/column_inputs.py @@ -0,0 +1,17 @@ +import dataclasses + +from ..constants.rebar_data import REBAR_SIZES, REBAR_STRENGTHS + + +@dataclasses.dataclass +class ColumnInputs: + w: float = 24 + h: float = 36 + bar_size: str = REBAR_SIZES[5] + bar_cover: float = 2 + bars_x: int = 6 + bars_y: int = 8 + fc: float = 8000 + fy: float = REBAR_STRENGTHS[1] + cover_to_center: bool = False + spiral_reinf: bool = False diff --git a/examples/conc_col_pmm/calc_document/dcr_calc_runner.py b/examples/conc_col_pmm/calc_document/dcr_calc_runner.py new file mode 100644 index 0000000..10213ca --- /dev/null +++ b/examples/conc_col_pmm/calc_document/dcr_calc_runner.py @@ -0,0 +1,32 @@ +from efficalc import FigureFromMatplotlib, Heading + +from ..calc_document.plotting import get_capacity, point_plotter +from ..col.axial_limits import AxialLimits +from ..col.column import Column +from ..pmm_search.ecc_search.get_dcr_ecc import get_dcr_ecc +from ..pmm_search.load_combo import LoadCombination + + +def calc_dcrs( + load_combos: list[LoadCombination], mesh, col: Column, axial_limits: AxialLimits +): + dcr_results = [] + for load in load_combos: + if load.show_in_report: # show the full calculations for this load case + Heading( + "DCR Calculation for Load Case P=" + + str(round(load.p, 1)) + + ", Mx=" + + str(round(load.mx, 1)) + + ", My=" + + str(round(load.my, 1)) + ) + capacity_pts = get_capacity.get_capacity(mesh, load) + # plot the PM diagram for this point + pm_figure = point_plotter.plot(capacity_pts, load, False) + FigureFromMatplotlib( + pm_figure, "PM interaction diagram for this load case. " + ) + + dcr_results.append(get_dcr_ecc(col, load, axial_limits)) + return dcr_results diff --git a/examples/conc_col_pmm/calc_document/document_wrapper.py b/examples/conc_col_pmm/calc_document/document_wrapper.py new file mode 100644 index 0000000..f141062 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/document_wrapper.py @@ -0,0 +1,42 @@ +from efficalc.report_builder import ReportBuilder + +from .calculation import calculation as col_input_calc + +# parameters: "override_inputs" is boolean and indicates whether the values +# provided in the next two arguments (column parameters, then load data) +# should override the default/user-input Input values + + +def run(override_inputs: bool, col_data: list, loads: list[list]): + if override_inputs: + new_inputs = {} + + input_to_name = { + "w": "w", + "h": "h", + "bar_size": "", + "bar_cover": " ", + "bars_x": " ", + "bars_y": " ", + "fc": r"f^{\prime}_c", + "fy": "f_y", + "cover_type": " ", + "transverse_type": " ", + } + for i in range(len(input_to_name)): + new_inputs[input_to_name["w"]] = col_data[0] + new_inputs[input_to_name["h"]] = col_data[1] + new_inputs[input_to_name["bar_size"]] = col_data[2] + new_inputs[input_to_name["bar_cover"]] = col_data[3] + new_inputs[input_to_name["bars_x"]] = col_data[4] + new_inputs[input_to_name["bars_y"]] = col_data[5] + new_inputs[input_to_name["fc"]] = col_data[6] + new_inputs[input_to_name["fy"]] = col_data[7] + new_inputs[input_to_name["cover_type"]] = col_data[8] + new_inputs[input_to_name["transverse_type"]] = col_data[9] + + builder = ReportBuilder(lambda: col_input_calc(default_loads=loads), new_inputs) + builder.view_report() + else: + builder = ReportBuilder(col_input_calc) + builder.view_report() diff --git a/examples/conc_col_pmm/calc_document/full_calc_document.py b/examples/conc_col_pmm/calc_document/full_calc_document.py new file mode 100644 index 0000000..23aaa92 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/full_calc_document.py @@ -0,0 +1,31 @@ +from ..calc_document.plotting import pure_mx_my_plotter +from ..col import assign_max_min +from ..col.col_canvas import draw_column_with_dimensions +from ..col.column import Column +from ..pmm_search.load_combo import LoadCombination +from .dcr_calc_runner import calc_dcrs +from .plotting.pmm_mesh import get_mesh +from .results_summary import results_summarizer + + +def calculation( + col: Column, + load_combos: list[LoadCombination], +): + # draw the column cross-section with dimensions and callouts + draw_column_with_dimensions.draw(col, "Section of Column") + + # calculate_axial_load_limits the max tension and compression to this column + axial_limits = assign_max_min.calculate_axial_load_limits(col) + + # Retrieve the quarter PMM mesh, which has points + # in the format (Mx, My, P). + _, _, _, mesh = get_mesh(col, 48, 18, axial_limits) + + # show the PM curves for bending purely about the x and y axes + pure_mx_my_plotter.plot(mesh) + + # calculate DCRs for all load cases + dcr_results = calc_dcrs(load_combos, mesh, col, axial_limits) + + results_summarizer(load_combos, dcr_results) diff --git a/examples/conc_col_pmm/calc_document/plotting/PMM.py b/examples/conc_col_pmm/calc_document/plotting/PMM.py new file mode 100644 index 0000000..ead8908 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/PMM.py @@ -0,0 +1,10 @@ +from ...pmm_search.load_combo import LoadCombination +import dataclasses + + +@dataclasses.dataclass +class PMM(): + X: list[float] + Y: list[float] + Z: list[float] + load_combos: list[LoadCombination] \ No newline at end of file diff --git a/examples/conc_col_pmm/calc_document/plotting/__init__.py b/examples/conc_col_pmm/calc_document/plotting/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/calc_document/plotting/get_capacity.py b/examples/conc_col_pmm/calc_document/plotting/get_capacity.py new file mode 100644 index 0000000..9b62a9b --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/get_capacity.py @@ -0,0 +1,54 @@ +import math + +from ...pmm_search.load_combo import LoadCombination + +""" +The function below interpolates between points on the PMM diagram to construct +the PM diagram for a given load point. It returns two lists which contain the +resultant moment and axial load for the points of the PM diagram at the lambda +for the given load point. +Parameters: "mesh" is the quarter PMM mesh consisting of (Mx, My, P) points and +"point" is the load point in the form (P, Mx, My). +""" + + +def get_capacity(mesh, point: LoadCombination): + pt_count = len(mesh) # the number of rows of points vertically + quarter = len(mesh[0]) - 1 # the number of angle spaces between points in a + # quadrant + + angle_space = (math.pi / 2) / quarter # the horizontal space between points + + lambda_transform = math.atan2(abs(point.my), abs(point.mx)) # the angle for + # the current point transformed to the range 0 to 90, where 0 must + # correspond to My=0 + + index = int(lambda_transform // angle_space) # the position of the mesh point in + # each row of "mesh" just below the lambda of the load point + + # this is how far beyond the chosen interval the load point lies + angle_extra = lambda_transform % angle_space + + # define factors that can be multiplied by the load values at two + # adjacent points to interpolate between those points + factors = [(angle_space - angle_extra) / angle_space, angle_extra / angle_space] + + phi_Mn = [0] * pt_count + phi_Pn = [0] * pt_count + + for i in range(pt_count): + # calculate the estimate of the biaxial moment and axial force capacity + # at the current point. The index for the angle of the current point must be + # limited to "quarter" to prevent it from exceeding the size of the + # mesh array, even if "index" is the max element due to rounding. + phi_Pn[i] = sum( + (mesh[i][min(index + j, quarter)][2] * factors[j] for j in range(2)) + ) + + for j in range(2): + # calculate the moment resultant for the given capacity point + moment = math.sqrt( + sum((mesh[i][min(index + j, quarter)][k] ** 2 for k in range(2))) + ) + phi_Mn[i] += moment * factors[j] + return [phi_Mn, phi_Pn] diff --git a/examples/conc_col_pmm/calc_document/plotting/get_pmm_data.py b/examples/conc_col_pmm/calc_document/plotting/get_pmm_data.py new file mode 100644 index 0000000..116528e --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/get_pmm_data.py @@ -0,0 +1,33 @@ +import numpy as np + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...pmm_search.load_combo import LoadCombination +from .pmm_mesh import get_mesh +from .PMM import PMM + +''' +This function takes inputs for a column and creates a +dataclass instance containing all the information for the +PMM diagram for the given column. That information is used +for plotting the PMM diagram. +''' + + +def get_pmm_data( + col: Column, + intervals: int, + load_spaces: int, + load_combos: list[LoadCombination], + axial_limits: AxialLimits, +): + # get the capacity point mesh for plotting the PMM surface. x, y, + # and z correspond to Mx, My, and P, respectively. "quarter_mesh" + # has just one quarter of the PMM mesh (including points aligned + # with the x and y axes) + x, y, z, _ = get_mesh(col, intervals, load_spaces, axial_limits) + X = np.array(x) + Y = np.array(y) + Z = np.array(z) + + return PMM(X=X, Y=Y, Z=Z, load_combos=load_combos) \ No newline at end of file diff --git a/examples/conc_col_pmm/calc_document/plotting/pmm_mesh.py b/examples/conc_col_pmm/calc_document/plotting/pmm_mesh.py new file mode 100644 index 0000000..a066025 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/pmm_mesh.py @@ -0,0 +1,106 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...pmm_search.load_search import bisect_load +from ...pmm_search.load_search.point_search_load import search + + +def get_mesh(col: Column, intervals, load_spaces, axial_limits: AxialLimits): + # "intervals" is the number of spaces in the angle of eccentricity, + # "load_spaces" is the number of vertical spaces in the PMM diagram + # Returns a mesh containing all the points of the PMM diagram, plus + # a quarter of that mesh + + # vectors containing the points to plot + x = [] + y = [] + z = [] + + # the bottom point must be added intervals+1 times because each level of + # the mesh needs intervals+1 points in order to form a closed surface + x.append([0] * (intervals + 1)) + y.append([0] * (intervals + 1)) + z.append([axial_limits.min_phi_pn] * (intervals + 1)) + + vert_space = axial_limits.load_span / load_spaces + c_guess = 0 # the starting guess for neutral axis + + quarter = math.floor(intervals / 4) # number of angle spaces in one quadrant + lambda_space = 2 * math.pi / intervals # the change in angle between points + theta_guesses = [-lambda_space * i for i in range(1, quarter)] # the guess for + # what theta should be at each point between the x and y axes, to be + # updated at each vertical level + count = 0 + + quarter_mesh = [[0] * (quarter + 1) for i in range(load_spaces + 2)] # a matrix + # to contain the output points, with height (outer dimension) spanning all + # the axial load interpolation points and width (inner dimension) spanning + # from lambda=0 to lambda=90 + quarter_mesh[0] = [(0, 0, axial_limits.min_phi_pn) for i in range(quarter + 1)] + quarter_mesh[load_spaces + 1] = [ + (0, 0, axial_limits.max_phi_pn) for i in range(quarter + 1) + ] + + for vert_count in range(1, load_spaces + 1): + count += 1 + c_guess += (col.w + col.h) / (2 * load_spaces) # increment the guess for + # the neutral axis depth + load_target = axial_limits.min_phi_pn + vert_space * vert_count + coords = [[0] * (intervals + 1) for i in range(3)] + + # get the output for theta=0 and apply it in 3 points + out = bisect_load.bisect(col, [0, load_target], [0, c_guess], axial_limits) + quarter_mesh[vert_count][0] = out[:3] + for pos in (0, quarter * 2, quarter * 4): + mult = [1 if (pos != quarter * 2) else -1, 1, 1] + for j in range(3): + coords[j][pos] = out[j] * mult[j] + + # get the output for theta=-90 and apply it in 2 points + out = bisect_load.bisect( + col, [math.pi / 2, load_target], [-math.pi / 2, c_guess], axial_limits + ) + quarter_mesh[vert_count][quarter] = out[:3] + for pos in (quarter, quarter * 3): + mult = [1, 1 if pos == quarter else -1, 1] + for j in range(3): + coords[j][pos] = out[j] * mult[j] + + # iterate over the remaining possible neutral axis angles between 0 + # and -90 and for each one, add its point to the four quadrants of the + # PMM diagram + lambda_target = lambda_space # the current target lambda + for i in range(1, quarter): + out = search( + col, + [lambda_target, load_target], + [theta_guesses[i - 1], c_guess], + axial_limits, + ) + quarter_mesh[vert_count][i] = out[:3] + theta_guesses[i - 1] = out[3] + + c_guess = out[4] + indices = [ + i, + quarter * 2 - i, + quarter * 2 + i, + intervals - i, + ] # the indices in + # the vector "coords" corresponding to the current point + for pos, index in enumerate(indices): + mult = [1 if (pos == 0 or pos == 3) else -1, 1 if pos < 2 else -1, 1] + for j in range(3): + coords[j][index] = out[j] * mult[j] + lambda_target += lambda_space + + x.append(coords[0]) + y.append(coords[1]) + z.append(coords[2]) + + x.append([0] * (intervals + 1)) + y.append([0] * (intervals + 1)) + z.append([axial_limits.max_phi_pn] * (intervals + 1)) + + return x, y, z, quarter_mesh diff --git a/examples/conc_col_pmm/calc_document/plotting/pmm_plotter_plotly.py b/examples/conc_col_pmm/calc_document/plotting/pmm_plotter_plotly.py new file mode 100644 index 0000000..68ac524 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/pmm_plotter_plotly.py @@ -0,0 +1,216 @@ +import numpy as np +import plotly.graph_objects as go + +""" +This function plots the factored load capacity diagram for the column +""" + + +def plot(pmm_data): + X = pmm_data.X + Y = pmm_data.Y + Z = pmm_data.Z + + axis_colors = {"x": "#FF0000", "y": "#00CC00", "z": "#0000FF"} + ax_labels = {"x": "${\phi}M_{nx}$", "y": "${\phi}M_{ny}$", "z": "${\phi}P_n$"} + + # factors for how far the axes should extend past the PMM surface + z_factor = 0.12 + xy_factor = 0.05 + + data = {} + # list for min, then max, then difference for each of X, Y, Z + min_max = [] + for dir in (X, Y, Z): + min_max.append([dir.min(), dir.max()]) + for i in range(3): + min_max[i].append(min_max[i][1] - min_max[i][0]) + min1, max1 = min_max[2][:2] + range1 = min_max[2][2] + data["z"] = {} + data["z"]["range"] = range1 * (1 + 2 * z_factor) + data["z"]["min"] = min1 - range1 * z_factor + data["z"]["max"] = max1 + range1 * z_factor + + # set the maximum aspect ratio (should be low so that the PMM + # surface won't be super long/short in any dimension + min_aspect = 2 + max_xy = max(min_max[0][2], min_max[1][2]) + for i, dir in enumerate(("x", "y")): + data[dir] = {} + length = min((1 + xy_factor) * max_xy, min_aspect * min_max[i][2]) + data[dir]["range"] = length + data[dir]["min"] = -length / 2 + data[dir]["max"] = length / 2 + + # creates an arrow which is used as the axis for whichever axis + # is input as the parameter + def get_arrow(axisname="x"): + ax_color = axis_colors[axisname] + scale = [[0, ax_color], [1, ax_color]] + # Create arrow body + body = go.Scatter3d( + marker=dict(size=1, color=ax_color), + line=dict(color=ax_color, width=5), + showlegend=False, # hide the legend + hoverinfo="skip", + ) + + head = go.Cone( + sizemode="raw", + sizeref=1, + autocolorscale=None, + colorscale=scale, + showscale=False, # disable additional colorscale for arrowheads + hoverinfo="skip", + ) + for ax, direction in zip(("x", "y", "z"), ("u", "v", "w")): + if ax == axisname: + body[ax] = data[ax]["min"], data[ax]["max"] + head[ax] = [data[ax]["max"]] + head[direction] = [0.05 * data[axisname]["range"]] + else: + body[ax] = 0, 0 + head[ax] = [0] + head[direction] = [0] + + return [body, head] + + def add_axis_arrows(fig): + for ax in ("x", "y", "z"): + for item in get_arrow(ax): + fig.add_trace(item) + + # returns a dictionary which forms the axis label for one of the axes + def get_annotation_for_ax(ax): + d = dict( + showarrow=False, + text=ax_labels[ax], + xanchor="left", + font=dict(color="#1f1f1f", size=18), + ) + for ax_ in ("x", "y", "z"): + if ax_ == ax: + d[ax_] = data[ax]["max"] - data[ax]["range"] * 0.025 + else: + d[ax_] = 0 + + if ax in {"x", "y"}: + d["xshift"] = 12 + else: + d["xshift"] = 2 + + return d + + def get_axis_names(): + return [get_annotation_for_ax(ax) for ax in ("x", "y", "z")] + + # returns the Plotly axes (should be empty since these default axes + # are not used) + def get_scene_axis(): + return dict( + title="", # remove axis label (x,y,z) + showbackground=False, + visible=False, + showticklabels=False, # hide numeric values of axes + showgrid=False, # Show box around plot + ) + + # plot the 3D PMM surface + fig = go.Figure( + layout=dict( + title={ + "text": "PMM Diagram", + "x": 0.5, + "y": 0.98, + "xanchor": "center", + "yanchor": "top", + "font": {"color": "#1f1f1f", "size": 18}, + }, + autosize=True, + width=550, + height=550, + margin=dict(l=5, r=5, b=5, t=5), + scene=dict( + xaxis=get_scene_axis(), + yaxis=get_scene_axis(), + zaxis=get_scene_axis(), + annotations=get_axis_names(), + aspectmode="cube", + ), + ), + ) + + add_axis_arrows(fig) + + # convert the PMM data to a numpy array for plotting + load_data = np.array([[ld.p, ld.mx, ld.my] for ld in pmm_data.load_combos]) + + # define a color (dark blue) for the load points + pt_color = "#002095" + + # Add points for the load cases. Note that because + # the load cases are passed in the form (P, Mx, My), + # the order must be changed for plotting + fig.add_trace( + go.Scatter3d( + mode="markers", + showlegend=False, # hide the legend + x=load_data[:, 1], + y=load_data[:, 2], + z=load_data[:, 0], + marker=dict( + color=pt_color, + size=4, + ), + ) + ) + + # define a color (orange) for the PMM surface + surface_color = "#ffbb0f" + surface_scale = [[0, surface_color], [1, surface_color]] + fig.add_trace( + go.Surface( + z=Z, + x=X, + y=Y, + opacity=0.5, + colorscale=surface_scale, + showscale=False, # Set to True to show colorscale + name="", + ) + ) + + # define a color (near-white) for the plotted mesh + line_color = "#f7f7f7" + + # plot a mesh on the PMM surface between the capacity points + line_size = 1.5 + for i in range(X.shape[0]): + fig.add_trace( + go.Scatter3d( + mode="lines", + line=dict(color=line_color, width=line_size), + showlegend=False, # hide the legend + x=X[i, :], + y=Y[i, :], + z=Z[i, :], + ) + ) + for i in range(X.shape[1]): + fig.add_trace( + go.Scatter3d( + mode="lines", + line=dict(color=line_color, width=line_size), + showlegend=False, # hide the legend + x=X[:-1, i], + y=Y[:-1, i], + z=Z[:-1, i], + ) + ) + fig.update_yaxes( + scaleanchor="x", + scaleratio=1, + ) + + return fig diff --git a/examples/conc_col_pmm/calc_document/plotting/point_plotter.py b/examples/conc_col_pmm/calc_document/plotting/point_plotter.py new file mode 100644 index 0000000..0988533 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/point_plotter.py @@ -0,0 +1,85 @@ +import math + +import matplotlib.pyplot as plt + +from ...pmm_search.load_combo import LoadCombination + +""" +The function "plot" creates and returns a matplotlib figure showing a PM diagram. + +Parameters: + +capacity_pts: a list of the capacity points on the PM curve, in the format +(Mxy, P). + +point: the load point if this plot is for a DCR calculation, or None if this plot +is simply to show the PM curve aligned with the Mx or My axis. The point is in +the format (P, Mx, My). + +only_Mx: boolean, and only used when point=None, True means this is for bending +in the x-direction. +""" + + +def plot(capacity_pts, point: LoadCombination | None, only_Mx): + [phi_Mn, phi_Pn] = capacity_pts + + if point: + # get the lamdba for the load point + pt_lambda = math.atan2(point.my, point.mx) # the angle for the current point + else: + # set the lambda to the desired axis + pt_lambda = 0 if only_Mx else math.pi / 2 + + pt_count = len(phi_Mn) # how many points are plotted + + fig, ax = plt.subplots() + ax.grid(True, which="both", zorder=0, linewidth=0.4) + # set the x-spine (see below for more info on `set_position`) + ax.spines["left"].set_position("zero") + + # turn off the right spine/ticks + ax.spines["right"].set_color("none") + ax.yaxis.tick_left() + + # set the y-spine + ax.spines["bottom"].set_position("zero") + + # turn off the top spine/ticks + ax.spines["top"].set_color("none") + ax.xaxis.tick_bottom() + + plot_angle = str(round(pt_lambda * 180 / math.pi)) + "$\degree$" + ax.set_title("P-M Interaction Diagram at $\lambda=$" + plot_angle) + + plt.plot(phi_Mn, phi_Pn, linewidth=2, zorder=1) + + ax.set_xlabel("${\phi}M_{nxy}$ (kip-ft)", fontsize=12, zorder=2) + ax.set_ylabel("${\phi}P_{n}$ (kip)", fontsize=12, zorder=2) + + ax.set_axisbelow(True) + + load_span = phi_Pn[-1] - phi_Pn[0] + label_offsets = (max(phi_Mn) * 0.008, load_span * 0.02) + # label the intersections with the y-axis + for i in (0, pt_count - 1): + pos = (phi_Mn[i], phi_Pn[i]) + label = str(round(phi_Pn[i], 1)) + plt.plot(pos[0], pos[1], marker="+", ms=12, mew=1.2, c="black", zorder=3) + plt.text(pos[0] + label_offsets[0], pos[1] + label_offsets[1], label, zorder=3) + + if point: + Muxy = math.sqrt(point.mx**2 + point.my**2) # the biaxial moment + + # plot and label the load point + pos = (Muxy, point.p) + moment_label = "($M_{uxy}=$" + str(round(Muxy, 1)) + " kip-ft, " + axial_label = "$P_u=$" + str(round(point.p, 1)) + " kip)" + label = moment_label + "\n" + axial_label + + plt.plot(pos[0], pos[1], marker="+", ms=12, mew=1.2, c="red", zorder=4) + plt.text(pos[0] + label_offsets[0], pos[1] + label_offsets[1], label, zorder=5) + + fig = ax.get_figure() + + return fig diff --git a/examples/conc_col_pmm/calc_document/plotting/pure_mx_my_plotter.py b/examples/conc_col_pmm/calc_document/plotting/pure_mx_my_plotter.py new file mode 100644 index 0000000..5a490cc --- /dev/null +++ b/examples/conc_col_pmm/calc_document/plotting/pure_mx_my_plotter.py @@ -0,0 +1,27 @@ +from efficalc import FigureFromMatplotlib, Heading + +from . import point_plotter + +""" +This function plots the intersection of the PMM surface with the vertical planes +Mx-P and My-P. +""" + + +def plot(mesh): + Heading("PM Diagrams for Pure Mx and My") + n = len(mesh) + m = len(mesh[0]) + capacity_pts = [ + [mesh[i][0][0] for i in range(n)], + [mesh[i][0][2] for i in range(n)], + ] + pm_figure = point_plotter.plot(capacity_pts, None, True) + FigureFromMatplotlib(pm_figure, "PM interaction diagram for pure Mx.") + + capacity_pts = [ + [mesh[i][m - 1][1] for i in range(n)], + [mesh[i][m - 1][2] for i in range(n)], + ] + pm_figure = point_plotter.plot(capacity_pts, None, False) + FigureFromMatplotlib(pm_figure, "PM interaction diagram for pure My.") diff --git a/examples/conc_col_pmm/calc_document/results_summary.py b/examples/conc_col_pmm/calc_document/results_summary.py new file mode 100644 index 0000000..595ea21 --- /dev/null +++ b/examples/conc_col_pmm/calc_document/results_summary.py @@ -0,0 +1,41 @@ +from efficalc import Calculation, Comparison, Heading, Table + +from ..pmm_search.load_combo import LoadCombination + +""" +Creates a table with the DCRs for all load cases. Also calculates the max +DCR, checks whether it is less than 1, and adds this to the result check. +""" + + +def results_summarizer(load_combos: list[LoadCombination], dcr_results): + Heading("Summary of Results") + data = [ + [ + ld.p, + ld.mx, + ld.my, + round(dcr_results[i], 2), + "O.K." if dcr_results[i] < 1 else "N.G.", + ] + for i, ld in enumerate(load_combos) + ] + + headers = ["Pu (kip)", "Mux (kip-ft)", "Muy (kip-ft)", "PM Vector DCR", "Passing?"] + Table(data, headers, "DCRs For All Load Cases", False, False) + + # calculate the max DCR and show + max_dcr = Calculation( + "DCR_{max}", + round(max(dcr_results), 2), + description="Maximum of all Load Case DCRs", + ) + Comparison( + max_dcr, + "<", + 1.0, + true_message="O.K.", + false_message="N.G.", + description="Max DCR check", + result_check=True, + ) diff --git a/examples/conc_col_pmm/calc_document/show_dcr_calc.py b/examples/conc_col_pmm/calc_document/show_dcr_calc.py new file mode 100644 index 0000000..813f3af --- /dev/null +++ b/examples/conc_col_pmm/calc_document/show_dcr_calc.py @@ -0,0 +1,68 @@ +from latexexpr_efficalc import Variable + +from efficalc import Calculation, Comparison, Heading, TextBlock, absolute, maximum + +from ..pmm_search.load_combo import LoadCombination +from .column_capacities import ColumnCapacities + + +def show(load: LoadCombination, capacity: ColumnCapacities): + + Heading(f"DCR Calculation - Load Case #{load.id}", 2) + + only_axial = capacity.Mx.result() == 0 and capacity.My.result() == 0 + if only_axial: + TextBlock( + "Since the load point is on the P axis, the DCR can be calculated by comparing the applied axial" + " load to the axial capacity calculated above:" + ) + p = Variable("P_u", load.p, "kip") + dcr = Calculation( + f"DCR_{{{load.id}}}", + absolute(p / capacity.P) if capacity.P.result() != 0 else 0, + ) + Comparison( + dcr, + "<", + 1.0, + true_message="O.K.", + false_message="N.G.", + description=f"Design check for load case #{load.id}", + ) + + else: + TextBlock( + "Compare the ratios of demand to capacity for Mx, My, and P to show that the calculated capacity" + " point is on the same PMM vector as the demand point. Note that the absolute value for the" + " moment DCRs is because the column has equal moment capacity in opposite directions by symmetry." + ) + + mux = Variable("M_{ux}", load.mx, "kip-ft") + dcr_mx = Calculation( + "DCR_{Mx}", absolute(mux / capacity.Mx) if capacity.Mx.result() != 0 else 0 + ) + + my = Variable("M_{uy}", load.my, "kip-ft") + dcr_my = Calculation( + "DCR_{My}", absolute(my / capacity.My) if capacity.My.result() != 0 else 0 + ) + + p = Variable("P_u", load.p, "kip") + dcr_p = Calculation( + "DCR_{P}", absolute(p / capacity.P) if capacity.P.result() != 0 else 0 + ) + + dcr = Calculation( + f"DCR_{{{load.id}}}", + maximum(dcr_mx, dcr_my, dcr_p), + "", + "The final DCR is:", + ) + Comparison( + dcr, + "<", + 1.0, + true_message="O.K.", + false_message="N.G.", + description=f"Design check for load case #{load.id}", + ) diff --git a/examples/conc_col_pmm/calc_document/try_axis_document.py b/examples/conc_col_pmm/calc_document/try_axis_document.py new file mode 100644 index 0000000..3d2c4ea --- /dev/null +++ b/examples/conc_col_pmm/calc_document/try_axis_document.py @@ -0,0 +1,626 @@ +import math + +from latexexpr_efficalc import Variable + +from efficalc import ( + PI, + Calculation, + ComparisonStatement, + Heading, + Symbolic, + Table, + TextBlock, + cos, + maximum, + minimum, + r_brackets, + sin, + tan, +) + +from ..col.axial_limits import AxialLimits +from ..col.col_canvas import draw_column_comp_zone, draw_column_with_triangle +from ..col.column import Column +from .column_capacities import ColumnCapacities + + +# Make sure that small numbers like 5.684e−14 get rounded to 0 +def _round(x: float): + return round(x, 5) + + +def try_axis_document( + col: Column, + axial_limits: AxialLimits, + theta_input=0, + c_input=10, +) -> ColumnCapacities: + + w = col.w_input + h = col.h_input + bar_area = col.rebar_area_input + fc = col.fc_input + fy = col.fy_input + E_s = col.steel_modulus_input + conc_epsilon = col.concrete_strain_input + + TextBlock( + "The neutral axis angle and depth below are iteratively determined to produce a capacity point aligning " + "exactly with the PMM vector of the applied load. \n" + ) + theta = Calculation( + "\\theta", _round(theta_input), "rad", description="Neutral axis angle" + ) + c = Calculation("c", _round(c_input), "in", description="Neutral axis depth") + + if c.get_value() == 0: + # This is dead code, this function doesn't get called for pure axial load cases, can probably remove + TextBlock( + "Because the neutral axis depth is zero, the column is in pure tension." + ) + return ColumnCapacities( + Variable("{\\phi}M_{nx}", 0, "kip-ft"), + Variable("{\\phi}M_{ny}", 0, "kip-ft"), + axial_limits.min_phi_pn_calculation, + ) + Heading("Forces in the Concrete", 2) + + if fc.get_value() <= 4000: + ComparisonStatement(2500, "<=", fc, "<=", 4000) + beta1 = Calculation( + "\\beta_1", + 0.85, + "", + "", + "ACI 318-19 Table 22.2.2.4.3(a)", + ) + elif fc.get_value() < 8000: + ComparisonStatement(4000, "<", fc, "<", 8000) + beta1 = Calculation( + "\\beta_1", + 0.85 - 0.05 * r_brackets(fc - 4000) / 1000, + "", + "", + "ACI 318-19 Table 22.2.2.4.3(b)", + ) + else: + ComparisonStatement(fc, ">=", 8000) + beta1 = Calculation( + "\\beta_1", + 0.65, + "", + "", + "ACI 318-19 Table 22.2.2.4.3(c)", + ) + a = Calculation( + "a", + beta1 * c, + "in", + "Depth of equivalent compression zone:", + "ACI 318-19 22.2.2.4.1", + ) + fc = Calculation( + "f^{\prime}_c", + fc / 1000, + "ksi", + description="Concrete strength converted to ksi:", + ) + epsilon = 1e-6 # acceptable error for considering the neutral axis to + # be vertical or horizontal + + intersects = [False] * 4 # whether line of concrete compression intersects the + # left, top, right, and bottom edges of the section + if theta.get_value() > -math.pi / 2 + epsilon: + left_y_temp = ( + col.half_h + - a.get_value() / math.cos(theta.get_value()) + - col.w * math.tan(theta.get_value()) + ) + intersects[0] = -col.half_h < left_y_temp < col.half_h + if intersects[0]: + left_y = Calculation( + "y_{\\mathrm{left}}", + h / 2 - a / cos(theta) - w * tan(theta), + "in", + " y coordinate of equivalent compression zone intersection with left edge:", + ) + + right_y = col.half_h - a.get_value() / math.cos(theta.get_value()) + intersects[2] = -col.half_h < right_y < col.half_h + if intersects[2]: + right_y = Calculation( + "y_{\\mathrm{right}}", + h / 2 - a / cos(theta), + "in", + " y coordinate of equivalent compression zone intersection with right edge:", + ) + + if theta.get_value() < -epsilon: + top_x_temp = col.half_w + a.get_value() / math.sin(theta.get_value()) + intersects[1] = -col.half_w < top_x_temp < col.half_w + if intersects[1]: + top_x = Calculation( + "x_{\\mathrm{top}}", + w / 2 + a / sin(theta), + "in", + " x coordinate of equivalent compression zone intersection with top edge:", + ) + + # if theta is -90, take the tan of 0 and get a change + # of zero, and if theta is -45, take the tan of 45 + # and get somewhat of an increase in x + bot_x_temp = ( + col.half_w + + a.get_value() / math.sin(theta.get_value()) + + col.h * math.tan(PI / 2 + theta.get_value()) + ) + intersects[3] = -col.half_w < bot_x_temp < col.half_w + if intersects[3]: + bot_x = Calculation( + "x_{\\mathrm{bottom}}", + w / 2 + a / sin(theta) + h * tan(PI / 2 + theta), + "in", + " x coordinate of equivalent compression zone intersection with bottom edge:", + ) + + # define accumulator variables + pn_tot = 0 + mnx_tot = 0 + mny_tot = 0 + + if not any(intersects): # the whole concrete section is in compression + TextBlock("The equivalent compression zone covers the whole concrete section. ") + pn_conc = Calculation("P_{\\mathrm{n, conc.}}", 0.85 * fc * w * h, "kips") + mnx_conc = Calculation( + "M_{\\mathrm{nx, conc.}}", + 0, + "kip-in", + ) + mny_conc = Calculation( + "M_{\\mathrm{ny, conc.}}", + 0, + "kip-in", + ) + else: + conc_area_num = 0 + + def add_axial_moment(pt_a, pt_b, pt_c): + nonlocal pn_tot, mnx_tot, mny_tot, conc_area_num + conc_area_num += 1 + Heading("Forces in Concrete Area " + str(conc_area_num), 3) + TextBlock( + "Below are the coordinates of the three points A, B, and C which define compression area" + " number " + str(conc_area_num) + "." + ) + pt_a = ( + Calculation("A_x", pt_a[0]), + Calculation("A_y", pt_a[1]), + ) + pt_b = ( + Calculation("B_x", pt_b[0]), + Calculation("B_y", pt_b[1]), + ) + pt_c = ( + Calculation("C_x", pt_c[0]), + Calculation("C_y", pt_c[1]), + ) + points = (pt_a, pt_b, pt_c) + points = [((pt[0]).get_value(), (pt[1]).get_value()) for pt in points] + draw_column_with_triangle.draw( + col, "Compression Area Outline", conc_area_num, points + ) + tri_area = Calculation( + "A_{\\mathrm{triangle}}", + 0.5 + * abs( + pt_a[0] * pt_b[1] + + pt_b[0] * pt_c[1] + + pt_c[0] * pt_a[1] + - pt_a[0] * pt_c[1] + - pt_b[0] * pt_a[1] + - pt_c[0] * pt_b[1] + ), + "in^2", + "Calculate the area of this triangular compression zone:", + ) + + centr_x = Calculation( + "x_{\\mathrm{centroid}}", + (pt_a[0] + pt_b[0] + pt_c[0]) / 3, + "in", + "x coordinate of the centroid of this zone:", + ) + centr_y = Calculation( + "y_{\\mathrm{centroid}}", + (pt_a[1] + pt_b[1] + pt_c[1]) / 3, + "in", + "y coordinate of the centroid of this zone:", + ) + return tri_area, centr_x, centr_y + + pt1 = (-w / 2, left_y) if intersects[0] else (top_x, h / 2) + pt2 = (bot_x, -h / 2) if intersects[3] else (w / 2, right_y) + points_block = [(col.w / 2, col.h / 2)] + if intersects[2]: + points_block.append((col.w / 2, right_y.get_value())) + else: + points_block.append((col.w / 2, -col.h / 2)) + points_block.append((bot_x.get_value(), -col.h / 2)) + if intersects[1]: + points_block.append((top_x.get_value(), col.h / 2)) + else: + points_block.append((-col.w / 2, left_y.get_value())) + points_block.append((-col.w / 2, col.h / 2)) + draw_column_comp_zone.draw( + col, "Equivalent compression zone outlined in red. ", points_block + ) + TextBlock( + "The equivalent stress block is now broken down into triangular areas and the forces are calculated for each." + ) + + (tri_area, centr_x, centr_y) = add_axial_moment( + pt1, pt2, (w / 2, h / 2) + ) # compression triangle to + # top right corner + pn_top_right = Calculation( + "P_{\\mathrm{n,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area, + "kips", + ) + mnx_top_right = Calculation( + "M_{\\mathrm{nx,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_y, + "kip-in", + ) + mny_top_right = Calculation( + "M_{\\mathrm{ny,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_x, + "kip-in", + ) + if intersects[0]: # there is a compression triangle to top left corner + (tri_area, centr_x, centr_y) = add_axial_moment( + (-w / 2, h / 2), (w / 2, h / 2), pt1 + ) + pn_top_left = Calculation( + "P_{\\mathrm{n,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area, + "kips", + ) + mnx_top_left = Calculation( + "M_{\\mathrm{nx,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_y, + "kip-in", + ) + mny_top_left = Calculation( + "M_{\\mathrm{ny,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_x, + "kip-in", + ) + if intersects[3]: # there is a compression triangle to bot right corner + (tri_area, centr_x, centr_y) = add_axial_moment( + (w / 2, -h / 2), (w / 2, h / 2), pt2 + ) + pn_bot_right = Calculation( + "P_{\\mathrm{n,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area, + "kips", + ) + mnx_bot_right = Calculation( + "M_{\\mathrm{nx,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_y, + "kip-in", + ) + mny_bot_right = Calculation( + "M_{\\mathrm{ny,\ Area\ " + str(conc_area_num) + "}}", + 0.85 * fc * tri_area * centr_x, + "kip-in", + ) + Heading("Total Forces in Concrete", 3) + if intersects[0] and intersects[3]: # take a sum for all 3 areas + pn_conc = Calculation( + "P_{\\mathrm{n, conc.}}", + pn_top_right + pn_top_left + pn_bot_right, + "kips", + ) + mnx_conc = Calculation( + "M_{\\mathrm{nx, conc.}}", + mnx_top_right + mnx_top_left + mnx_bot_right, + "kip-in", + ) + mny_conc = Calculation( + "M_{\\mathrm{ny, conc.}}", + mny_top_right + mny_top_left + mny_bot_right, + "kip-in", + ) + elif intersects[0]: + pn_conc = Calculation( + "P_{\\mathrm{n, conc.}}", pn_top_right + pn_top_left, "kips" + ) + mnx_conc = Calculation( + "M_{\\mathrm{nx, conc.}}", + mnx_top_right + mnx_top_left, + "kip-in", + ) + mny_conc = Calculation( + "M_{\\mathrm{ny, conc.}}", + mny_top_right + mny_top_left, + "kip-in", + ) + elif intersects[3]: + pn_conc = Calculation( + "P_{\\mathrm{n, conc.}}", pn_top_right + pn_bot_right, "kips" + ) + mnx_conc = Calculation( + "M_{\\mathrm{nx, conc.}}", + mnx_top_right + mnx_bot_right, + "kip-in", + ) + mny_conc = Calculation( + "M_{\\mathrm{ny, conc.}}", + mny_top_right + mny_bot_right, + "kip-in", + ) + else: + pn_conc = Calculation("P_{\\mathrm{n, conc.}}", pn_top_right, "kips") + mnx_conc = Calculation( + "M_{\\mathrm{nx, conc.}}", + mnx_top_right, + "kip-in", + ) + mny_conc = Calculation( + "M_{\\mathrm{ny, conc.}}", + mny_top_right, + "kip-in", + ) + + Heading("Equations for Rebar Axial and Moment Calculations", 2) + + TextBlock( + "Each bar is at coordinates (x,y) relative to the column centroid. For example, the top right bar" + " is located at the coordinates below:" + ) + right_bar_x = col.half_w - col.edge_to_bar_center # x coordinate of bars on the + # right edge + top_bar_y = col.half_h - col.edge_to_bar_center # y coordinate of bars on the + # top edge + + x = Calculation("x_{\mathrm{bar}}", right_bar_x) + y = Calculation("y_{\mathrm{bar}}", top_bar_y) + + d_bar = Symbolic( + "d_{\mathrm{bar}}", + r_brackets(w / 2 - x) * cos(theta + PI / 2) + + r_brackets(h / 2 - y) * sin(theta + PI / 2), + "Effective depth:", + ) + strain_sym = Symbolic( + "\epsilon_{\mathrm{bar}}", conc_epsilon * (d_bar - c) / c, "Strain:" + ) + stress_sym = Symbolic( + "\sigma_{\mathrm{bar}}", + minimum(fy, maximum(-fy, conc_epsilon * strain_sym)), + "Stress:", + ) + Symbolic( + "\sigma_{\mathrm{bar}}", + stress_sym + 0.85 * fc, + "If the bar is in the equivalent compression zone, add the concrete stress to avoid double-counting:", + ) + TextBlock( + "Note that the sign of the following expressions is reversed because positive stress in the rebar " + "is defined as tension while positive axial moment in the column is defined as compression." + ) + Symbolic( + "P_{\mathrm{bar}}", + -bar_area * stress_sym * y, + "Contribution to moment about the " "x axis:", + ) + Symbolic( + "M_{\mathrm{x, bar}}", + -bar_area * stress_sym * y, + "Contribution to moment about " "the x axis:", + ) + Symbolic( + "M_{\mathrm{y, bar}}", + -bar_area * stress_sym * x, + "Contribution to moment about " "the y axis:", + ) + # define new accumulators + pn = 0 + mnx = 0 + mny = 0 + + bar_num = 0 + + def add_bar(coords): + nonlocal pn, mnx, mny, bar_num + bar_num += 1 + bar_calc = [bar_num] + bar_calc.extend([round(val, 2) for val in coords]) + # "offset" is the distance from the center of the bar to the line + # passing through the top right corner of the section and parallel to + # the neutral axis + offset = (col.w / 2 - coords[0]) * math.cos(theta.get_value() + math.pi / 2) + ( + col.h / 2 - coords[1] + ) * math.sin(theta.get_value() + math.pi / 2) + bar_calc.append(round(offset, 2)) + strain = conc_epsilon.get_value() * (offset - c.get_value()) / c.get_value() + bar_calc.append(round(strain, 4)) + stress = min( + fy.get_value(), + max( + -fy.get_value(), + conc_epsilon.get_value() + * E_s.get_value() + * (offset - c.get_value()) + / c.get_value(), + ), + ) + bar_calc.append(round(stress, 2)) + in_comp_zone = offset < a.get_value() + if in_comp_zone: + stress += 0.85 * fc.get_value() + bar_calc.append(round(0.85 * fc.get_value(), 2)) + else: + bar_calc.append(0) + # since negative strain and negative stress are defined as + # compression for rebar but compression is positive in the conc. + # the sign of everything needs to be changed + pn -= bar_area.get_value() * stress + bar_calc.append(round(-bar_area.get_value() * stress, 1)) + mnx -= bar_area.get_value() * stress * coords[1] + # the +0 is to avoid rounding to -0 + bar_calc.append(round(-bar_area.get_value() * stress * coords[1], 1) + 0) + mny -= bar_area.get_value() * stress * coords[0] + bar_calc.append(round(-bar_area.get_value() * stress * coords[0], 1) + 0) + + rebar_matrix.append(bar_calc) + + y = col.y_start + # iterate over the bars along the left and right lines + # (this includes corner bars) + bar_count = 0 + + headers = [ + "Bar Number", + "X Coord. (in)", + "Y Coord. (in)", + "Effective Depth d (in)", + "Strain (unitless)", + "Stress (ksi)", + "Stress Correction for Displaced Concrete (ksi)", + "Axial Force (kips)", + "Contribution to Mx (kip-in)", + "Contribution to My (kip-in)", + ] + rebar_matrix = [] + for i in range(col.bars_y): + for x in (-right_bar_x, right_bar_x): + coords = [x, y] + bar_count += 1 + add_bar(coords) + y += col.y_space + + x = col.x_start + for i in range(col.bars_x - 2): + # iterate over the bars along the top and bottom lines, and add the + # force for each one + for y in (-top_bar_y, top_bar_y): + coords = [x, y] + add_bar(coords) + x += col.x_space + Heading("Tabulated Rebar Results", 2) + + Table(rebar_matrix, headers) + + Heading("Force Totals", 2) + pn_steel = Calculation("P_{\\mathrm{n, steel}}", pn, "kips") + mnx_steel = Calculation( + "M_{\\mathrm{nx, steel}}", + _round(mnx), + "kip-in", + ) + mny_steel = Calculation( + "M_{\\mathrm{ny, steel}}", + _round(mny), + "kip-in", + ) + pn = Calculation("P_{\\mathrm{n, tot}}", pn_conc + pn_steel, "kips") + mnx = Calculation( + "M_{\\mathrm{nx, tot}}", + mnx_conc + mnx_steel, + "kip-in", + ) + mny = Calculation( + "M_{\\mathrm{ny, tot}}", + mny_conc + mny_steel, + "kip-in", + ) + + Heading("Capacity Calculation", 2) + + TextBlock("The extreme tension reinforcement is centered at these coordinates:") + coords = [0, 0] + coords[0] = Calculation("x_{\\mathrm{bar}}", -right_bar_x, "in") + coords[1] = Calculation("y_{\\mathrm{bar}}", -top_bar_y, "in") + offset = Calculation( + "d_t", + r_brackets(w / 2 - coords[0]) * cos(theta + PI / 2) + + r_brackets(h / 2 - coords[1]) * sin(theta + PI / 2), + "in", + ) + max_strain = Calculation( + "\\epsilon_y", conc_epsilon * r_brackets(offset - c) / c, "" + ) + yield_strain = Calculation("\\epsilon_{ty}", fy / E_s, "") + if max_strain.get_value() <= yield_strain.get_value(): + ComparisonStatement(max_strain, "<=", yield_strain) + strain_level = 0 + elif ( + yield_strain.get_value() + < max_strain.get_value() + < yield_strain.get_value() + 0.003 + ): + ComparisonStatement(yield_strain, "<", max_strain, "<", yield_strain + 0.003) + strain_level = 1 + else: + ComparisonStatement(max_strain, ">=", yield_strain + 0.003) + strain_level = 2 + if col.spiral_reinf: + if strain_level == 0: + phi = Calculation( + "\\phi", + 0.75, + "", + "Failure is compression-controlled and transverse reinforcement is spiral:", + "ACI 318-19 Table 21.2.2(a)", + ) + elif strain_level == 1: + phi = Calculation( + "\\phi", + 0.75 + 0.15 * r_brackets(max_strain - yield_strain) / 0.003, + "", + "Failure is transition and transverse reinforcement is spiral:", + "ACI 318-19 Table 21.2.2(c)", + ) + else: + phi = Calculation( + "\\phi", + 0.90, + "", + "Failure is tension-controlled and transverse reinforcement is spiral:", + "ACI 318-19 Table 21.2.2(e)", + ) + else: + if strain_level == 0: + phi = Calculation( + "\\phi", + 0.65, + "", + "Failure is compression-controlled and transverse reinforcement is tied:", + "ACI 318-19 Table 21.2.2(b)", + ) + elif strain_level == 1: + phi = Calculation( + "\\phi", + 0.65 + 0.25 * r_brackets(max_strain - yield_strain) / 0.003, + "", + "Failure is transition and transverse reinforcement is tied:", + "ACI 318-19 Table 21.2.2(d)", + ) + else: + phi = Calculation( + "\\phi", + 0.90, + "", + "Failure is tension-controlled and transverse reinforcement is tied:", + "ACI 318-19 Table 21.2.2(f)", + ) + TextBlock("Factored axial and moment capacities:") + phi_pn = Calculation( + "{\\phi}P_n", minimum(phi * pn, axial_limits.max_phi_pn_calculation), "kips" + ) + phi_mnx = Calculation("{\\phi}M_{nx}", phi * mnx / 12, "kip-ft") + phi_mny = Calculation("{\\phi}M_{ny}", phi * mny / 12, "kip-ft") + + return ColumnCapacities(phi_mnx, phi_mny, phi_pn) diff --git a/examples/conc_col_pmm/calc_report_examples/calc_report_example1.pdf b/examples/conc_col_pmm/calc_report_examples/calc_report_example1.pdf new file mode 100644 index 0000000..2443ee3 Binary files /dev/null and b/examples/conc_col_pmm/calc_report_examples/calc_report_example1.pdf differ diff --git a/examples/conc_col_pmm/calc_report_examples/calc_report_example2.pdf b/examples/conc_col_pmm/calc_report_examples/calc_report_example2.pdf new file mode 100644 index 0000000..0b2957b Binary files /dev/null and b/examples/conc_col_pmm/calc_report_examples/calc_report_example2.pdf differ diff --git a/examples/conc_col_pmm/col/__init__.py b/examples/conc_col_pmm/col/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/col/assign_max_min.py b/examples/conc_col_pmm/col/assign_max_min.py new file mode 100644 index 0000000..b2f00cc --- /dev/null +++ b/examples/conc_col_pmm/col/assign_max_min.py @@ -0,0 +1,95 @@ +from efficalc import Calculation, Heading, TextBlock, r_brackets + +from ..struct_analysis import try_axis +from .axial_limits import AxialLimits +from .column import Column + + +def calculate_axial_load_limits(col: Column) -> AxialLimits: + + w = col.w_input + h = col.h_input + bar_area = col.rebar_area_input + bars_x = col.bars_x_input + bars_y = col.bars_y_input + fc = col.fc_input + fy = col.fy_input + + Heading("Axial Capacity Calculations") + + steel_area = Calculation( + "A_{st}", + bar_area * r_brackets(2 * bars_x + 2 * bars_y - 4), + "in^2", + "Total area of longitudinal reinforcement:", + ) + tot_area = Calculation("A_g", w * h, "in^2", "Gross section area") + Heading("Compressive Capacity", 2) + max_pn = Calculation( + "P_0", + 0.85 * fc / 1000 * r_brackets(tot_area - steel_area) + fy * steel_area, + "kips", + "", + "ACI 318-19 22.4.2.2", + ) + if col.spiral_reinf: + TextBlock("Because the transverse reinforcement is spiral:") + max_pn_limit = Calculation( + "P_{\mathrm{n,max}}", + 0.85 * max_pn, + "kips", + "Maximum axial strength", + "ACI 318-19 22.4.2.1(b)", + ) + phi = Calculation( + "\\phi", + 0.75, + "", + "", + "ACI 318-19 Table 21.2.2(a)", + ) + else: + TextBlock("Because the transverse reinforcement is tied:") + max_pn_limit = Calculation( + "P_{\mathrm{n,max}}", + 0.80 * max_pn, + "kips", + "", + "ACI 318-19 22.4.2.1(a)", + ) + phi = Calculation( + "\\phi", + 0.65, + "", + "", + "ACI 318-19 Table 21.2.2(b)", + ) + + max_phi_pn = Calculation("{\\phi}P_{\mathrm{n,max}}", phi * max_pn_limit, "kips") + + Heading("Tensile Capacity", 2) + min_pn = Calculation( + "P_{\mathrm{nt,max}}", + -1 * fy * steel_area, + "kips", + "", + "ACI 318-19 22.4.3.1", + ) + + phi = Calculation( + "\\phi", + try_axis.PHI_FLEXURE, + "", + "Because failure is tension-controlled:", + "ACI 318-19 21.2.2(e)", + ) + min_phi_pn = Calculation("{\\phi}P_{\mathrm{nt,max}}", phi * min_pn, "kips") + + return AxialLimits( + max_pn.result(), + max_phi_pn.result(), + max_phi_pn, + min_pn.result(), + min_phi_pn.result(), + min_phi_pn, + ) diff --git a/examples/conc_col_pmm/col/axial_limits.py b/examples/conc_col_pmm/col/axial_limits.py new file mode 100644 index 0000000..8249c64 --- /dev/null +++ b/examples/conc_col_pmm/col/axial_limits.py @@ -0,0 +1,19 @@ +import dataclasses + +from efficalc import Calculation + + +@dataclasses.dataclass +class AxialLimits: + max_pn: float + max_phi_pn: float + max_phi_pn_calculation: Calculation + min_pn: float + min_phi_pn: float + min_phi_pn_calculation: Calculation + + # difference between the maximum and minimum allowable loads, + # to be used for normalizing error + @property + def load_span(self) -> float: + return self.max_phi_pn - self.min_phi_pn diff --git a/examples/conc_col_pmm/col/col_canvas/__init__.py b/examples/conc_col_pmm/col/col_canvas/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/col/col_canvas/draw_column_comp_zone.py b/examples/conc_col_pmm/col/col_canvas/draw_column_comp_zone.py new file mode 100644 index 0000000..c6f1341 --- /dev/null +++ b/examples/conc_col_pmm/col/col_canvas/draw_column_comp_zone.py @@ -0,0 +1,76 @@ +import math + +from efficalc.canvas import Canvas, Dimension, Polyline, Rectangle, Text + +from ..column import Column +from .draw_plain_column import draw as draw_base + + +def draw(col: Column, caption_input: str, points) -> Canvas: + # number is the zone number with which to label this triangle, points are the three + # corners of the triangle as a tuple of 2-element tuples + canvas = draw_base(col) + canvas.caption = caption_input + + points = [(pt[0] + col.w / 2, col.h / 2 - pt[1]) for pt in points] + n = len(points) + centr = [sum((points[i][j] for i in range(n))) / n for j in range(2)] + points.append(points[0]) + + canvas.add(Polyline(points, stroke_width=0.2, stroke="#c80c00")) + + canvas.add( + Rectangle( + centr[0] - 0.3, + centr[1] - 1, + 10, + 1.4, + 0.2, + 0.2, + stroke_width=0.12, + stroke="#c80c00", + fill="white", + ) + ) + canvas.add( + Text( + "Compression Zone", + centr[0], + centr[1], + font_size=1.2, + ) + ) + + # margin around the section + m = 8 + scale_factor = 0.37817187 * math.log((col.w + col.h) / 2) + 0.03808133 + # add dimensions + common_dim_styles = { + "unit": '"', + "gap": 0.15, + "stroke_width": 0.08 * scale_factor, + "text_size": 1.2, + } + for pt in points[:-1]: + if pt[0] == 0 and 0 < pt[1] < col.h: + # this point is on the left edge, add dimension + canvas.add(Dimension(0, 0, 0, pt[1], offset=-0.5 * m, **common_dim_styles)) + if pt[0] == col.w and 0 < pt[1] < col.h: + # this point is on the right edge, add dimension + canvas.add( + Dimension(col.w, 0, col.w, pt[1], offset=0.5 * m, **common_dim_styles) + ) + if pt[1] == 0 and 0 < pt[0] < col.w: + # this point is on the top edge, add dimension + canvas.add( + Dimension(col.w, 0, pt[0], 0, offset=0.5 * m, **common_dim_styles) + ) + if pt[1] == col.h and 0 < pt[0] < col.w: + # this point is on the bottom edge, add dimension + canvas.add( + Dimension( + col.w, col.h, pt[0], col.h, offset=-0.5 * m, **common_dim_styles + ) + ) + + return canvas diff --git a/examples/conc_col_pmm/col/col_canvas/draw_column_with_dimensions.py b/examples/conc_col_pmm/col/col_canvas/draw_column_with_dimensions.py new file mode 100644 index 0000000..2801b56 --- /dev/null +++ b/examples/conc_col_pmm/col/col_canvas/draw_column_with_dimensions.py @@ -0,0 +1,117 @@ +import math + +from efficalc.canvas import ArrowMarker, Canvas, Dimension, Leader, Line, Text + +from ...constants.rebar_data import rebar_diameter +from ..column import Column +from .draw_plain_column import draw as draw_base + + +def draw(col: Column, caption_input: str, unit: str = '"') -> Canvas: + canvas = draw_base(col) + canvas.caption = caption_input + scale_factor = 0.37817187 * math.log((col.w + col.h) / 2) + 0.03808133 + + # margin around the section + m = 8 + + # x and y axes + canvas.add( + Line( + col.w / 2, + col.h / 2, + col.w + m / 4, + col.h / 2, + stroke="black", + stroke_width=0.08, + marker_end=ArrowMarker(), + ) + ) + canvas.add( + Line( + col.w / 2, + col.h / 2, + col.w / 2, + -m / 4, + stroke="black", + stroke_width=0.08, + marker_end=ArrowMarker(), + ) + ) + + # define reinforcement properties + bars_x = col.bars_x + bars_y = col.bars_y + long_bar_radius = rebar_diameter(col.bar_size) / 2 + + # define "cover" to be the cover to center of bar (regardless of whether the user specified the + # cover to be clear or to center) + cover = col.bar_cover if col.cover_to_center else col.bar_cover + long_bar_radius + + canvas.add(Text("x", col.w + m / 4 + 0.25, col.h / 2, font_size=1)) + canvas.add(Text("y", col.w / 2, -m / 4 - 0.25, font_size=1)) + + # add dimensions + common_dim_styles = { + "unit": unit, + "gap": 0.15, + "stroke_width": 0.08 * scale_factor, + "text_size": 1.2, + } + canvas.add(Dimension(0, 0, col.w, 0, offset=0.5 * m, **common_dim_styles)) + canvas.add(Dimension(col.w, 0, col.w, col.h, offset=0.5 * m, **common_dim_styles)) + canvas.add( + Dimension( + 0, + cover, + col.bar_cover, + cover, + offset=0.25 * m + cover, + **common_dim_styles, + ) + ) + + # add leaders + common_leader_styles = { + "marker": ArrowMarker(), + "landing_len": 1, + "stroke_width": 0.08 * scale_factor, + "text_size": 1.2, + } + + x_bar_starting_x = cover + x_bar_spacing = (col.w - 2 * cover) / (bars_x - 1) + x_bar_y_bot = col.h - cover + + y_bar_starting_y = x_bar_starting_x + y_bar_spacing = (col.h - 2 * cover) / (bars_y - 1) + y_bar_x_left = x_bar_starting_x + + bottom_bar_x = x_bar_starting_x + (col.bars_x - 2) * x_bar_spacing + bottom_bar_y = x_bar_y_bot + canvas.add( + Leader( + bottom_bar_x + long_bar_radius * 0.85, + bottom_bar_y + long_bar_radius * 0.85, + bottom_bar_x + cover + m / 8, + col.h + m / 8, + f"({col.bars_x}){col.bar_size} x-direction, E.S.", + **common_leader_styles, + ) + ) + left_bar_x = y_bar_x_left + long_bar_radius * 0.85 + left_bar_y = ( + y_bar_starting_y + (col.bars_y - 2) * y_bar_spacing + long_bar_radius * 0.85 + ) + canvas.add( + Leader( + left_bar_x, + left_bar_y, + cover + 2 * long_bar_radius + x_bar_spacing, + col.h + m / 5, + f"({col.bars_y}){col.bar_size} y-direction, N.S.", + **common_leader_styles, + ) + ) + + return canvas diff --git a/examples/conc_col_pmm/col/col_canvas/draw_column_with_triangle.py b/examples/conc_col_pmm/col/col_canvas/draw_column_with_triangle.py new file mode 100644 index 0000000..9473a31 --- /dev/null +++ b/examples/conc_col_pmm/col/col_canvas/draw_column_with_triangle.py @@ -0,0 +1,100 @@ +import math + +from efficalc.canvas import Canvas, CircleMarker, Dimension, Polyline, Rectangle, Text + +from ...struct_analysis import triangles +from ..column import Column +from .draw_plain_column import draw as draw_base + + +def draw(col: Column, caption_input: str, number, points) -> Canvas: + # number is the zone number with which to label this triangle, points are the three + # corners of the triangle as a tuple of 2-element tuples + canvas = draw_base(col) + canvas.caption = caption_input + + points = [(pt[0] + col.w / 2, col.h / 2 - pt[1]) for pt in points] + points.append(points[0]) + + marker = CircleMarker() + canvas.add( + Polyline( + points, + stroke_width=0.2, + stroke="#c80c00", + marker_mid=marker, + marker_end=marker, + ) + ) + + centr = triangles.triangle_centroid(*points[:3]) + canvas.add( + Rectangle( + centr[0] - 0.3, + centr[1] - 1, + 3.8, + 1.4, + 0.2, + 0.2, + stroke_width=0.12, + stroke="#c80c00", + fill="white", + ) + ) + canvas.add( + Text( + "Area " + str(number), + centr[0], + centr[1], + font_size=1.2, + ) + ) + + pt_labels = ["Point A", "Point B", "Point C"] + for i in range(3): + offsets = ( + [0.2] * 2 + if (points[i][0] > 1e-6 and points[i][1] < col.h - 1e-6) + else [-3, -1.1] + ) + canvas.add( + Text( + pt_labels[i], + points[i][0] + offsets[0], + points[i][1] - offsets[1], + font_size=1, + ) + ) + + # margin around the section + m = 8 + scale_factor = 0.37817187 * math.log((col.w + col.h) / 2) + 0.03808133 + # add dimensions + common_dim_styles = { + "unit": '"', + "gap": 0.15, + "stroke_width": 0.08 * scale_factor, + "text_size": 1.2, + } + for pt in points: + if pt[0] == 0 and 0 < pt[1] < col.h: + # this point is on the left edge, add dimension + canvas.add(Dimension(0, 0, 0, pt[1], offset=-0.5 * m, **common_dim_styles)) + if pt[0] == col.w and 0 < pt[1] < col.h: + # this point is on the right edge, add dimension + canvas.add( + Dimension(col.w, 0, col.w, pt[1], offset=0.5 * m, **common_dim_styles) + ) + if pt[1] == 0 and 0 < pt[0] < col.w: + # this point is on the top edge, add dimension + canvas.add( + Dimension(col.w, 0, pt[0], 0, offset=0.5 * m, **common_dim_styles) + ) + if pt[1] == col.h and 0 < pt[0] < col.w: + # this point is on the bottom edge, add dimension + canvas.add( + Dimension( + col.w, col.h, pt[0], col.h, offset=-0.5 * m, **common_dim_styles + ) + ) + return canvas diff --git a/examples/conc_col_pmm/col/col_canvas/draw_plain_column.py b/examples/conc_col_pmm/col/col_canvas/draw_plain_column.py new file mode 100644 index 0000000..74cd71b --- /dev/null +++ b/examples/conc_col_pmm/col/col_canvas/draw_plain_column.py @@ -0,0 +1,106 @@ +from efficalc.canvas import Canvas, Circle, Rectangle + +from ...constants.rebar_data import rebar_diameter + + +def draw(col) -> Canvas: + # margin around the section + m = 8 + + # define reinforcement properties + bars_x = col.bars_x + bars_y = col.bars_y + long_bar_radius = rebar_diameter(col.bar_size) / 2 + stirrup_diameter = rebar_diameter(col.shear_bar_size) + stirrup_bend_radius = 3 * stirrup_diameter + + # define "cover" to be the cover to center of bar (regardless of whether the user specified the + # cover to be clear or to center) + cover = col.bar_cover if col.cover_to_center else col.bar_cover + long_bar_radius + + # set up the canvas + canvas = Canvas( + col.w + 2 * m, + col.h + 2 * m, + min_xy=(-m, -m), + scale=20, + default_element_stroke_width=0, + ) + + # Draw the beam outline + column_outline = Rectangle(0, 0, col.w, col.h, fill="rgb(0 0 0 / 30%)") + canvas.add(column_outline) + + # add transverse reinforcement + stirrups = Rectangle( + cover - long_bar_radius - stirrup_diameter / 2, + cover - long_bar_radius - stirrup_diameter / 2, + col.w - 2 * cover + 2 * long_bar_radius + stirrup_diameter, + col.h - 2 * cover + 2 * long_bar_radius + stirrup_diameter, + rx=stirrup_bend_radius, + ry=stirrup_bend_radius, + stroke_width=stirrup_diameter, + stroke="#004aad", + ) + + canvas.add(stirrups) + + # add longitudinal reinforcement + x_bar_starting_x = cover + x_bar_spacing = (col.w - 2 * cover) / (bars_x - 1) + x_bar_y_bot = col.h - cover + x_bar_y_top = x_bar_starting_x + + y_bar_starting_y = x_bar_starting_x + y_bar_spacing = (col.h - 2 * cover) / (bars_y - 1) + y_bar_x_left = x_bar_starting_x + y_bar_x_right = col.w - cover + + for i in range(bars_x): + # bottom bar + canvas.add( + Circle( + x_bar_starting_x + i * x_bar_spacing, + x_bar_y_bot, + long_bar_radius, + fill="black", + ) + ) + + # top bar + canvas.add( + Circle( + x_bar_starting_x + i * x_bar_spacing, + x_bar_y_top, + long_bar_radius, + fill="black", + ) + ) + + last_y_bar = bars_y - 1 + for i in range(bars_y): + if i == 0 or i == last_y_bar: + # corner bars are drawn as x bar + pass + + # left bar + canvas.add( + Circle( + y_bar_x_left, + y_bar_starting_y + i * y_bar_spacing, + long_bar_radius, + fill="black", + ) + ) + + # right bar + canvas.add( + Circle( + y_bar_x_right, + y_bar_starting_y + i * y_bar_spacing, + long_bar_radius, + fill="black", + ) + ) + + return canvas diff --git a/examples/conc_col_pmm/col/column.py b/examples/conc_col_pmm/col/column.py new file mode 100644 index 0000000..1f55879 --- /dev/null +++ b/examples/conc_col_pmm/col/column.py @@ -0,0 +1,104 @@ +from efficalc import Calculation, Input + +from ..constants.rebar_data import STEEL_E, rebar_area, rebar_diameter + +""" +The class Column takes the following parameters: + w = section width (x-dir) (in.) + h = section height (y-dir) (in.) + bar_size = rebar number (Imperial) + bar_cover = concrete cover to the longitudinal rebar (in.) + bars_x = number of bars on the top/bottom edge + bars_y = number of bars on the left/right edge + fc = concrete f'c (psi) + fy = steel yield strength (ksi) + cover_to_center = boolean: concrete cover is to the center of the rebar, + false means it is clear cover (to edge of bars) + spiral_reinf = boolean: the shear reinforcement is spiral + rebar_area = area of one longitudinal bar (in^2) + steel_modulus = steel modulus of elasticity (ksi) + concrete_strain = max concrete strain at f'c +""" + + +class Column: + # a constant for now + shear_bar_size = "#4" + + def __init__( + self, + w_input: Input, + h_input: Input, + bar_size_input: Input, + bar_cover_input: Input, + bars_x_input: Input, + bars_y_input: Input, + fc_input: Input, + fy_input: Input, + cover_to_center: bool, + spiral_reinf: bool, + rebar_area_input: Calculation, + steel_modulus_input: Calculation, + concrete_strain_input: Calculation, + ): + + self.w_input = w_input + self.h_input = h_input + self.bar_size_input = bar_size_input + self.bar_cover_input = bar_cover_input + self.bars_x_input = bars_x_input + self.bars_y_input = bars_y_input + self.fc_input = fc_input + self.fy_input = fy_input + self.rebar_area_input = rebar_area_input + self.steel_modulus_input = steel_modulus_input + self.concrete_strain_input = concrete_strain_input + + self.cover_to_center = cover_to_center + self.spiral_reinf = spiral_reinf + + self.w = w_input.get_value() + self.h = h_input.get_value() + self.bar_size = bar_size_input.get_value() + self.bar_cover = bar_cover_input.get_value() + self.bars_x = bars_x_input.get_value() + self.bars_y = bars_y_input.get_value() + self.fc = fc_input.get_value() + self.fy = fy_input.get_value() + + self.half_w = self.w / 2 + self.half_h = self.h / 2 + self.area = self.w * self.h + + # coordinates of three corner points + self.bot_right = (self.half_w, -self.half_h) + self.top_right = (self.half_w, self.half_h) + self.top_left = (-self.half_w, self.half_h) + + self.bar_area = rebar_area(self.bar_size) + self.bar_dia = rebar_diameter(self.bar_size) + if cover_to_center: # the concrete cover is to center of bars + self.edge_to_bar_center = self.bar_cover + else: + self.edge_to_bar_center = self.bar_cover + self.bar_dia / 2 + + # the center-to-center spacing of the top/bottom bars + self.x_space = (self.w - 2 * self.edge_to_bar_center) / (self.bars_x - 1) + + # the center-to-center spacing of the left/right bars + self.y_space = (self.h - 2 * self.edge_to_bar_center) / (self.bars_y - 1) + + # the y-coordinate of the bottom left corner bar + self.y_start = -self.half_h + self.edge_to_bar_center + + # the x-coordinate of the second-from-left bar on the bottom edge + self.x_start = -self.half_w + self.edge_to_bar_center + self.x_space + + self.beta1 = max( + 0.65, min(0.85, 0.85 - 0.05 / 1000 * (self.fc - 4000)) + ) # the ratio + # of a/c for this column + self.steel_yield = self.fy / STEEL_E # strain in steel at yielding + + # safety factor for compression-controlled column + self.PHI_COMP = 0.75 if spiral_reinf else 0.65 diff --git a/examples/conc_col_pmm/constants/__init__.py b/examples/conc_col_pmm/constants/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/constants/concrete_data.py b/examples/conc_col_pmm/constants/concrete_data.py new file mode 100644 index 0000000..1b7b275 --- /dev/null +++ b/examples/conc_col_pmm/constants/concrete_data.py @@ -0,0 +1 @@ +MAX_CONCRETE_STRAIN = 0.003 diff --git a/examples/conc_col_pmm/constants/rebar_data.py b/examples/conc_col_pmm/constants/rebar_data.py new file mode 100644 index 0000000..ea3f09a --- /dev/null +++ b/examples/conc_col_pmm/constants/rebar_data.py @@ -0,0 +1,17 @@ +from typing import Literal + +REBAR_SIZES = ["#3", "#4", "#5", "#6", "#7", "#8", "#9", "#10", "#11", "#14", "#18"] +REBAR_DIAMETERS = [0.38, 0.50, 0.63, 0.75, 0.88, 1.00, 1.13, 1.27, 1.41, 1.69, 2.26] +REBAR_AREAS = [0.11, 0.20, 0.31, 0.44, 0.60, 0.79, 1.00, 1.27, 1.56, 2.25, 4.00] +STEEL_E = 29000 # steel modulus of elasticity in ksi +REBAR_STRENGTHS = [40, 60, 80] + +BarSize = Literal["#3", "#4", "#5", "#6", "#7", "#8", "#9", "#10", "#11", "#14", "#18"] + + +def rebar_area(size: BarSize): + return REBAR_AREAS[REBAR_SIZES.index(size)] + + +def rebar_diameter(size: BarSize): + return REBAR_DIAMETERS[REBAR_SIZES.index(size)] diff --git a/examples/conc_col_pmm/images/ColumnNA-cropped.svg b/examples/conc_col_pmm/images/ColumnNA-cropped.svg new file mode 100644 index 0000000..546c07b --- /dev/null +++ b/examples/conc_col_pmm/images/ColumnNA-cropped.svg @@ -0,0 +1,104 @@ + + + + + + + + + + + + + + + + + + + + x + + + + + + y + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + θ + N.A. + + + \ No newline at end of file diff --git a/examples/conc_col_pmm/images/minimizer-cropped.svg b/examples/conc_col_pmm/images/minimizer-cropped.svg new file mode 100644 index 0000000..c26db80 --- /dev/null +++ b/examples/conc_col_pmm/images/minimizer-cropped.svg @@ -0,0 +1,115 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + f=0 + g=0 + f + + + + + + g + + + + + + Next iteration point + + + + + + Current iteration point + + + + + + + + \ No newline at end of file diff --git a/examples/conc_col_pmm/pmm_search/__init__.py b/examples/conc_col_pmm/pmm_search/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/__init__.py b/examples/conc_col_pmm/pmm_search/ecc_search/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/change_ecc.py b/examples/conc_col_pmm/pmm_search/ecc_search/change_ecc.py new file mode 100644 index 0000000..c0cc836 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/ecc_search/change_ecc.py @@ -0,0 +1,38 @@ +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...struct_analysis import try_axis +from .get_error_ecc import get_error + +delta = 1e-8 # small change to be used for finite differences + + +# this function is the same as "change" except that it finds a direction for +# eccentricity rather than for c +def change(col: Column, guess, target, output, axial_limits: AxialLimits): + error = get_error(output, target) + # A small positive value "delta" is added to both inputs in order to test + # the effect on the results of "try_axis". It may have to be negative for + # theta to avoid exceeding 0. + delta0 = delta if guess[0] < -delta else -delta + + output2 = try_axis.try_axis(col, guess[0] + delta0, guess[1], axial_limits) + a = (output2[0] - output[0]) / delta0 + c = (output2[1] - output[1]) / delta0 + + output2 = try_axis.try_axis(col, guess[0], guess[1] + delta, axial_limits) + b = (output2[0] - output[0]) / delta + d = (output2[1] - output[1]) / delta + + e = target[0] - output[0] + f = target[1] - output[1] + + change = [0] * 2 + det = a * d - b * c + + # avoid divide by zero + if det != 0: + # set the planned changes in theta and c to try to reach the point at + # which lambda and load are both their target values + change[0] = (d * e - b * f) / det + change[1] = (a * f - c * e) / det + return change, error diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/get_dcr_ecc.py b/examples/conc_col_pmm/pmm_search/ecc_search/get_dcr_ecc.py new file mode 100644 index 0000000..57cbb77 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/ecc_search/get_dcr_ecc.py @@ -0,0 +1,75 @@ +import math + +from latexexpr_efficalc import Variable + +from ...calc_document.column_capacities import ColumnCapacities +from ...calc_document.show_dcr_calc import show as show_dcr_calc +from ...calc_document.try_axis_document import try_axis_document +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ..load_combo import LoadCombination +from . import point_search_ecc + + +# accepts as arguments the column and the load point and returns +# the dcr for this particular load point. +def get_dcr_ecc(col: Column, load: LoadCombination, axial_limits: AxialLimits): + target_M = math.sqrt(load.mx**2 + load.my**2) + target_ecc = math.atan2(load.p, target_M) + + tol = 0.01 * math.pi / 180 + # Load case is in pure tension + if target_ecc > math.pi / 2 - tol: + if load.show_in_report: + capacities = ColumnCapacities( + Variable("{\\phi}M_{nx}", 0, "kip-ft"), + Variable("{\\phi}M_{ny}", 0, "kip-ft"), + axial_limits.max_phi_pn_calculation, + ) + show_dcr_calc(load, capacities) + return load.p / axial_limits.max_phi_pn + + # Load case is in pure compression + if target_ecc < -math.pi / 2 + tol: + if load.show_in_report: + capacities = ColumnCapacities( + Variable("{\\phi}M_{nx}", 0, "kip-ft"), + Variable("{\\phi}M_{ny}", 0, "kip-ft"), + axial_limits.min_phi_pn_calculation, + ) + show_dcr_calc(load, capacities) + return load.p / axial_limits.min_phi_pn + + # Load has a bending moment component + target_lambda = math.atan2(abs(load.my), abs(load.mx)) + target = (target_lambda, target_ecc) + + # the best guess of theta is -lambda + guess = [ + -target_lambda, + (col.w + col.h) / 2, + ] + + # find the point on the PMM diagram that is on the same vector as the + # applied load + Mx, My, P, final_guess = point_search_ecc.search(col, target, guess, axial_limits) + + dcrs = [float("inf")] * 3 + dcr_candidates = [] + if P != 0: + dcrs[0] = load.p / P + dcr_candidates.append(load.p / P) + if Mx != 0: + dcrs[1] = abs(load.mx / Mx) + dcr_candidates.append(abs(load.mx / Mx)) + if My != 0: + dcrs[2] = abs(load.my / My) + dcr_candidates.append(abs(load.my / My)) + dcr = max(dcr_candidates) if len(dcr_candidates) > 0 else 0 + if load.show_in_report: + capacities = try_axis_document( + col, axial_limits, final_guess[0], final_guess[1] + ) + + show_dcr_calc(load, capacities) + return dcr diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/get_error_ecc.py b/examples/conc_col_pmm/pmm_search/ecc_search/get_error_ecc.py new file mode 100644 index 0000000..94a448a --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/ecc_search/get_error_ecc.py @@ -0,0 +1,11 @@ +import math + + +def get_error(output, target): + # find the difference between both outputs and their target values + # print("in get error") + # print(output) + # print(target) + lambda_diff = output[0] - target[0] + ecc_diff = output[1] - target[1] + return math.sqrt(lambda_diff**2 + ecc_diff**2) diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/limit_comp_ecc.py b/examples/conc_col_pmm/pmm_search/ecc_search/limit_comp_ecc.py new file mode 100644 index 0000000..185d77e --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/ecc_search/limit_comp_ecc.py @@ -0,0 +1,28 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...struct_analysis import try_axis +from . import get_error_ecc + + +def limit_comp(col: Column, guess, target, axial_limits: AxialLimits): + guess[0] = min(0, max(guess[0], -math.pi / 2)) + guess[1] = max(1e-6, guess[1]) + # guess[1]=max(c_lims[0],max(c_lims[1],guess[1])) + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + lim_factor = 0.999 + while output[3] > lim_factor * col.PHI_COMP * axial_limits.max_pn: + # the current phi_pn (without the 0.8) is at or almost at its + # maximum value, which means the column is probably in full + # compression, which must be avoided or derivatives will be zero + guess[1] /= 2 + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + while output[3] < lim_factor * axial_limits.min_phi_pn: + # the current phi_pn is at or almost at its minimum value, so c must + # b increased + guess[1] *= 10 + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + # update the distance from the target point + error = get_error_ecc.get_error(output, target) + return output, error diff --git a/examples/conc_col_pmm/pmm_search/ecc_search/point_search_ecc.py b/examples/conc_col_pmm/pmm_search/ecc_search/point_search_ecc.py new file mode 100644 index 0000000..5024636 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/ecc_search/point_search_ecc.py @@ -0,0 +1,66 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from . import change_ecc, limit_comp_ecc + +""" +Searches for a point on the PMM diagram with a certain lambda and eccentricity. +Parameters: the column object, the target (as a tuple of lambda, then axial load), and the +initial guess (as a tuple of theta, then neutral axis depth) +Returns: Mx, My, P, and the final guess as a two-value list +*the angle of eccentricity guess must be between 0 and -pi/2 and the target lambda +must be between 0 and pi/2 +""" + + +# the purpose of this function is to search for a point with a particular lambda and eccentricity +def search(col: Column, target, guess, axial_limits: AxialLimits): + tol = 0.001 # error accepted as the actual point + + # get output for the initial guess point + output, error = limit_comp_ecc.limit_comp(col, guess, target, axial_limits) + + count = 1 # count of iterations (calls to "try_axis") + count_lim = 100 # iteration limit + + while error > tol and count < count_lim: + # get a descent direction for the current point + direction, error = change_ecc.change(col, guess, target, output, axial_limits) + + # since finding the direction requires two "try_axis" calls + count += 2 + + # the factor to be applied to the change direction + factor = 1 + + # try the guess point and decrease "factor" until the error is less + # than the error from the last point + error2 = error + 1 + while error2 > error and factor > 0.01: + guess2 = [guess[i] + factor * direction[i] for i in range(2)] + + output, error2 = limit_comp_ecc.limit_comp( + col, guess2, target, axial_limits + ) + # if "limit_comp" resulted in a change in the guess of c, update + # the change factor to save calls to "try_axis" in the next + # iteration. Since "guess2" was passed to "limit_comp", it is + # already updated + if guess2[1] != guess[1] + factor * direction[1] and direction[1] != 0: + factor = (guess2[1] - guess[1]) / direction[1] + + count += 1 + factor *= 0.6 + + guess = guess2 + error = error2 + + # return the forces at the final trial point + Mx = output[4] * math.cos(output[0]) + My = output[4] * math.sin(output[0]) + + # it is possible that the point will be on the top plateau, so the + # axial force must be limited + P = min(axial_limits.max_phi_pn, output[3]) + return Mx, My, P, guess diff --git a/examples/conc_col_pmm/pmm_search/load_combo.py b/examples/conc_col_pmm/pmm_search/load_combo.py new file mode 100644 index 0000000..e5360cb --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_combo.py @@ -0,0 +1,15 @@ +import dataclasses + + +@dataclasses.dataclass +class LoadCombination: + id: int + p: float + mx: float + my: float + show_in_report: bool + + +def is_yes(show: str | float): + trimmed_yes = f"{show}".strip().lower() + return trimmed_yes == "yes" or trimmed_yes == "y" diff --git a/examples/conc_col_pmm/pmm_search/load_search/__init__.py b/examples/conc_col_pmm/pmm_search/load_search/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/pmm_search/load_search/bisect_load.py b/examples/conc_col_pmm/pmm_search/load_search/bisect_load.py new file mode 100644 index 0000000..27560d7 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/bisect_load.py @@ -0,0 +1,89 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from .limit_comp_load import limit_comp +from .starting_pts import starting_pts + + +def bisect(col: Column, target, guess, axial_limits: AxialLimits): + # This function returns a point on the PMM diagram (Mx, My, and P) plus the + # two inputs that produced that point (theta and c) as a tuple. The point + # returned is intended to match the values in "target," which is a + # list containing the target lambda and target factored axial load in that + # order. "col" is the column being analyzed and "guess" contains the + # starting guess for c. In this function, c is the only input varied, + # and theta is held constant. This is because this function is intended + # for points on the PMM diagram where it is known that lambda=-theta, + # which occur on the positive and negative x and y axes. + + tol = 0.005 # normalized error accepted as the actual point + + depth = col.w * math.sin(target[0]) + col.h * math.cos(target[0]) # the distance + # between the section corners perpendicular to the neutral axis + + # set the two initial guess points + pts = starting_pts(col, guess, depth, target, axial_limits) + + error = float("inf") # normalized distance of the current point from the + # target + + guess = 0 # the guess for c at each iteration + change = 0 # the change in c between iterations + best_error = 10 # record for normalized error, initialize to a large number + best = [] # the point encountered so far with smallest normalized error + + # debug: + points = [] + points.append(pts[0]) + points.append(pts[1]) + + counter = 0 + while error > tol and counter < 50: + slope = (pts[1][3] - pts[0][3]) / (pts[1][1] - pts[0][1]) + + # calculate the distance of this output from its target + dist = target[1] - pts[1][3] + if slope != 0: # check slope to avoid dividing by zero + # move by the amount estimated to get to zero, minus a small + # reduction for stability + change = dist / slope + else: + change = 0 + + error = float("inf") + factor = 1 # factor for reducing the change between iterations + + # if the error increases after the first guess, the change should be + # reduced to try to improve the guess + while error > best_error and factor >= 0.1 and counter < 50: + # calculate the next guess based on the current point + guess = pts[1][1] + change * factor + factor /= 10 + # set the floor to avoid guesses that are zero or negative and to + # avoid setting the guess equal to the current point + guess_floor = 1e-6 if pts[1][1] != 1e-6 else 0.1 + guess = max(guess_floor, guess) + + # get the output for this guess, reducing c if the compression + # is too high to + output, error = limit_comp( + col, [-target[0]] + [guess], target, axial_limits + ) + counter += 1 + + # update the current and previous point + pts[0] = pts[1] + pts[1] = [-target[0]] + [guess] + [output[0]] + [output[3]] + points.append(pts[1]) + + # update the lowest error and the best point if this guess is a record + if error < best_error: + best_error = error + # the factored moments and axial force, plus the best guess inputs + Mx = output[4] * math.cos(output[0]) + My = output[4] * math.sin(output[0]) + P = output[3] + best = (Mx, My, P, pts[1][0], pts[1][1]) + + return best diff --git a/examples/conc_col_pmm/pmm_search/load_search/change_load.py b/examples/conc_col_pmm/pmm_search/load_search/change_load.py new file mode 100644 index 0000000..054e29b --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/change_load.py @@ -0,0 +1,36 @@ +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...struct_analysis import try_axis +from .get_error_load import get_error + +delta = 1e-8 # small change to be used for finite differences + + +def change(col: Column, guess, target, output, axial_limits: AxialLimits): + error = get_error(output, target, axial_limits.load_span) + # A small positive value "delta" is added to both inputs in order to test + # the effect on the results of "try_axis". It may have to be negative for + # theta to avoid exceeding 0. + delta0 = delta if guess[0] < -delta else -delta + + output2 = try_axis.try_axis(col, guess[0] + delta0, guess[1], axial_limits) + a = (output2[0] - output[0]) / delta0 + c = (output2[3] - output[3]) / delta0 + + output2 = try_axis.try_axis(col, guess[0], guess[1] + delta, axial_limits) + b = (output2[0] - output[0]) / delta + d = (output2[3] - output[3]) / delta + + e = target[0] - output[0] + f = target[1] - output[3] + + change = [0] * 2 + det = a * d - b * c + + # avoid divide by zero + if det != 0: + # set the planned changes in theta and c to try to reach the point at + # which lambda and load are both their target values + change[0] = (d * e - b * f) / det + change[1] = (a * f - c * e) / det + return change, error diff --git a/examples/conc_col_pmm/pmm_search/load_search/get_error_load.py b/examples/conc_col_pmm/pmm_search/load_search/get_error_load.py new file mode 100644 index 0000000..8b11c81 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/get_error_load.py @@ -0,0 +1,8 @@ +import math + + +def get_error(output, target, load_span): + # find the difference between both outputs and their target values + lambda_diff = output[0] - target[0] + load_diff = output[3] - target[1] + return math.sqrt((lambda_diff / (math.pi / 2)) ** 2 + (load_diff / load_span) ** 2) diff --git a/examples/conc_col_pmm/pmm_search/load_search/limit_comp_load.py b/examples/conc_col_pmm/pmm_search/load_search/limit_comp_load.py new file mode 100644 index 0000000..0516a70 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/limit_comp_load.py @@ -0,0 +1,28 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from ...struct_analysis import try_axis +from ..load_search.get_error_load import get_error + + +def limit_comp(col: Column, guess, target, axial_limits: AxialLimits): + guess[0] = min(0, max(guess[0], -math.pi / 2)) + guess[1] = max(1e-6, guess[1]) + # guess[1]=max(c_lims[0],max(c_lims[1],guess[1])) + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + lim_factor = 0.999 + while output[3] > lim_factor * col.PHI_COMP * axial_limits.max_pn: + # the current phi_pn (without the 0.8) is at or almost at its + # maximum value, which means the column is probably in full + # compression, which must be avoided or derivatives will be zero + guess[1] /= 2 + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + while output[3] < lim_factor * axial_limits.min_phi_pn: + # the current phi_pn is at or almost at its minimum value, so c must + # b increased + guess[1] *= 10 + output = try_axis.try_axis(col, guess[0], guess[1], axial_limits) + # update the distance from the target point + error = get_error(output, target, axial_limits.load_span) + return output, error diff --git a/examples/conc_col_pmm/pmm_search/load_search/point_search_load.py b/examples/conc_col_pmm/pmm_search/load_search/point_search_load.py new file mode 100644 index 0000000..3a199c7 --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/point_search_load.py @@ -0,0 +1,62 @@ +import math + +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from .change_load import change +from .limit_comp_load import limit_comp + +""" +Searches for a point on the PMM diagram with a certain axial load capacity and lambda. +Parameters: the column object, the target (as a tuple of lambda, then axial load), and the +initial guess (as a tuple of theta, then neutral axis depth) +Returns: Mx, My, P, and the two final iteration values for the final guess +*the angle of eccentricity guess must be between 0 and -pi/2 and the target lambda +must be between 0 and pi/2 +""" + + +def search(col: Column, target, guess, axial_limits: AxialLimits): + tol = 0.001 # error accepted as the actual point + + # get output for the initial guess point + output, error = limit_comp(col, guess, target, axial_limits) + + count = 1 # count of iterations (calls to "try_axis") + count_lim = 100 # iteration limit + + while error > tol and count < count_lim: + # get a descent direction for the current point + direction, error = change(col, guess, target, output, axial_limits) + + # since finding the direction requires two "try_axis" calls + count += 2 + + # the factor to be applied to the change direction + factor = 1 + + # try the guess point and decrease "factor" until the error is less + # than the error from the last point + error2 = error + 1 + while error2 > error and factor > 0.01: + guess2 = [guess[i] + factor * direction[i] for i in range(2)] + + output, error2 = limit_comp(col, guess2, target, axial_limits) + + # if "limit_comp" resulted in a change in the guess of c, update + # the change factor to save calls to "try_axis" in the next + # iteration. Since "guess2" was passed to "limit_comp", it is + # already updated + if guess2[1] != guess[1] + factor * direction[1] and direction[1] != 0: + factor = (guess2[1] - guess[1]) / direction[1] + + count += 1 + factor *= 0.6 + + guess = guess2 + error = error2 + + # return the forces at the final trial point + Mx = output[4] * math.cos(output[0]) + My = output[4] * math.sin(output[0]) + P = output[3] + return Mx, My, P, guess[0], guess[1] diff --git a/examples/conc_col_pmm/pmm_search/load_search/starting_pts.py b/examples/conc_col_pmm/pmm_search/load_search/starting_pts.py new file mode 100644 index 0000000..3f56d0d --- /dev/null +++ b/examples/conc_col_pmm/pmm_search/load_search/starting_pts.py @@ -0,0 +1,29 @@ +from ...col.axial_limits import AxialLimits +from ...col.column import Column +from .limit_comp_load import limit_comp + +reduction = 0.005 # the fraction of the total estimated span of both inputs +# that should be added/subtracted to the starting guess points + + +# Returns a list of two points which are the starting guesses for a gradient +# descent problem, where the two points are close together, centered on the +# supplied guess "guess," and different in both their theta and c values. +# "depth" is an estimate of the maximum c, and load_only is boolean, where +# True indicates that only the load should be varied +def starting_pts(col: Column, guess, depth, target, axial_limits: AxialLimits): + # calculate the starting differences in theta and c between two points + c_change = reduction * depth + change_factors = (-1, 1) # factors used to decrease parameters for A and to + # increase parameters for B + + pts = [guess.copy() for i in range(2)] + + for i, factor in enumerate(change_factors): + pts[i][1] += c_change * factor + guess[1] = max(1e-6, guess[1]) + + # calculate and store the load output for the current guess point + output, error = limit_comp(col, pts[i], target, axial_limits) + pts[i].extend([output[0], output[3]]) + return pts diff --git a/examples/conc_col_pmm/requirements.txt b/examples/conc_col_pmm/requirements.txt new file mode 100644 index 0000000..7c5f19b --- /dev/null +++ b/examples/conc_col_pmm/requirements.txt @@ -0,0 +1,5 @@ +pytest==8.0.2 +plotly==5.24.1 +efficalc +matplotlib +pandas \ No newline at end of file diff --git a/examples/conc_col_pmm/struct_analysis/__init__.py b/examples/conc_col_pmm/struct_analysis/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/struct_analysis/triangles.py b/examples/conc_col_pmm/struct_analysis/triangles.py new file mode 100644 index 0000000..54195d0 --- /dev/null +++ b/examples/conc_col_pmm/struct_analysis/triangles.py @@ -0,0 +1,13 @@ +def triangle_area(pt_a, pt_b, pt_c): + return 0.5 * abs( + pt_a[0] * pt_b[1] + + pt_b[0] * pt_c[1] + + pt_c[0] * pt_a[1] + - pt_a[0] * pt_c[1] + - pt_b[0] * pt_a[1] + - pt_c[0] * pt_b[1] + ) + + +def triangle_centroid(pt_a, pt_b, pt_c): + return ((pt_a[0] + pt_b[0] + pt_c[0]) / 3, (pt_a[1] + pt_b[1] + pt_c[1]) / 3) diff --git a/examples/conc_col_pmm/struct_analysis/try_axis.py b/examples/conc_col_pmm/struct_analysis/try_axis.py new file mode 100644 index 0000000..d784af2 --- /dev/null +++ b/examples/conc_col_pmm/struct_analysis/try_axis.py @@ -0,0 +1,157 @@ +import math +import random + +from ..col.axial_limits import AxialLimits +from ..col.column import Column +from .triangles import * + +PHI_FLEXURE = 0.9 # safety factor for flexure-controlled column +COMP_FACTOR = 0.8 # additional reduction factor for axial compression + + +def try_axis(col: Column, theta, c, axial_limits: AxialLimits): + # this function returns the lambda, eccentricity, pn, phi_pn, and phi_mn + # from particular neutral axis angle and neutral axis depth c. The neutral + # axis angle must be between -90 degrees and 0 degrees, inclusive, and the + # neutral axis depth must be greater than or equal to 0 + + if c == 0: + return 0, 0, axial_limits.min_pn, axial_limits.min_phi_pn, 0 + a = col.beta1 * c + epsilon = 1e-11 # acceptable error for considering the neutral axis to + # be vertical or horizontal + red_fc = 0.85 * col.fc / 1000 # reduced f'c for concrete compression limit, ksi + pn = 0 # total axial force, kips (positive is compression) + mn = [0, 0] # Mnx, then Mny, kip-in (positive is compression to the right/top) + + intersects = [False] * 4 # whether line of concrete compression intersects the + # left, top, right, and bottom edges of the section + if theta > -math.pi / 2 + epsilon: + left_y = col.half_h - a / math.cos(theta) - col.w * math.tan(theta) + intersects[0] = -col.half_h < left_y < col.half_h + + right_y = col.half_h - a / math.cos(theta) + intersects[2] = -col.half_h < right_y < col.half_h + + if theta < -epsilon: + top_x = col.half_w + a / math.sin(theta) + intersects[1] = -col.half_w < top_x < col.half_w + + # if theta is -90, take the tan of 0 and get a change + # of zero, and if theta is -45, take the tan of 45 + # and get somewhat of an increase in x + bot_x = col.half_w + a / math.sin(theta) + col.h * math.tan(math.pi / 2 + theta) + intersects[3] = -col.half_w < bot_x < col.half_w + + if not any(intersects): # the whole concrete section is in compression + pn += red_fc * col.area + else: + + def add_axial_moment(pt_a, pt_b, pt_c): + nonlocal pn + tri_area = triangle_area(pt_a, pt_b, pt_c) + pn += tri_area * red_fc + centr = triangle_centroid(pt_a, pt_b, pt_c) + for i in range(2): + mn[i] += tri_area * red_fc * centr[1 - i] + + pt1 = (-col.half_w, left_y) if intersects[0] else (top_x, col.half_h) + pt2 = (bot_x, -col.half_h) if intersects[3] else (col.half_w, right_y) + + add_axial_moment(pt1, pt2, col.top_right) # compression triangle to + # top right corner + if intersects[0]: # there is a compression triangle to top left corner + add_axial_moment(col.top_left, col.top_right, pt1) + if intersects[3]: # there is a compression triangle to bot right corner + add_axial_moment(col.bot_right, col.top_right, pt2) + + strain_per_in = -col.concrete_strain_input.get_value() / c + steel_max_strain = 0 # value to keep record of greatest steel tensile strain + + # calculate a normal vector rotated 90 degrees from the neutral axis angle + normal = (math.cos(theta + math.pi / 2), math.sin(theta + math.pi / 2)) + + def add_bar(coords): + nonlocal steel_max_strain + nonlocal pn + # "offset" is the distance from the center of the bar to the line + # passing through the top right corner of the section and parallel to + # the neutral axis + offset = normal[0] * (col.half_w - coords[0]) + normal[1] * ( + col.half_h - coords[1] + ) + strain = (c - offset) * strain_per_in + steel_max_strain = max(steel_max_strain, strain) + + stress = strain * col.steel_modulus_input.get_value() + stress = max(-col.fy, min(col.fy, stress)) + if a > offset: + # this means this bar is within the compression range, + # so subtract the stress in the concrete, the stress + # will be negative in this case, so add to it + stress += red_fc + # since negative strain and negative stress are defined as + # compression for rebar but compression is positive in the conc. + # the sign of everything needs to be changed + force = col.bar_area * stress + pn -= force + for i in range(2): + mn[i] -= force * coords[1 - i] + + right_bar_x = col.half_w - col.edge_to_bar_center # x coordinate of bars on the + # right edge + y = col.y_start + # iterate over the bars along the left and right lines + # (this includes corner bars) + for i in range(col.bars_y): + for x in (-right_bar_x, right_bar_x): + coords = (x, y) + add_bar(coords) + y += col.y_space + + top_bar_y = col.half_h - col.edge_to_bar_center # y coordinate of bars on the + # top edge + x = col.x_start + for i in range(col.bars_x - 2): + # iterate over the bars along the top and bottom lines, and add the + # force for each one + for y in (-top_bar_y, top_bar_y): + coords = (x, y) + add_bar(coords) + x += col.x_space + + lambda1 = math.atan2(mn[1], mn[0]) # atan2 takes y,x, and we want mny at top + + # if Mx and My are both zero, the angle isn't defined, so return a random + # number to help avoid divide by zero errors + if not any(mn): + lambda1 = random.uniform(0, math.pi / 2) + mn_xy = math.sqrt(mn[0] ** 2 + mn[1] ** 2) / 12 # the moment resultant in kip-ft + + # calculate the factor of safety depending on the maximum steel strain + phi = min( + PHI_FLEXURE, + max( + col.PHI_COMP, + col.PHI_COMP + + (PHI_FLEXURE - col.PHI_COMP) + * (steel_max_strain - col.steel_yield) + / col.concrete_strain_input.get_value(), + ), + ) + # the following values ignore the limit on phi_pn to help convergence near that value + # (so there will still be derivatives above that value) + phi_pn_not_limited = phi * pn + phi_mn_xy = phi * mn_xy + + # the eccentricity does include the limit on phi_pn, and it is reported as + # angle from the Mx-My axis for numerical stability. + phi_pn = min(phi_pn_not_limited, axial_limits.max_phi_pn) + ecc = math.atan2( + phi_pn, phi_mn_xy + ) # the eccentricity as an angle above the M-M plane + + # returning the angle of eccentricity, the eccentricity as an angle from Mx-My plane, the nominal + # axial capacity in kip, the factored axial capacity in kip, and the + # factored moment capacity in kip-ft + return lambda1, ecc, pn, phi_pn_not_limited, phi_mn_xy diff --git a/examples/conc_col_pmm/tests/__init__.py b/examples/conc_col_pmm/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/tests/conftest.py b/examples/conc_col_pmm/tests/conftest.py new file mode 100644 index 0000000..46ea9e3 --- /dev/null +++ b/examples/conc_col_pmm/tests/conftest.py @@ -0,0 +1,93 @@ +import pytest + +from efficalc import Calculation, Input + +from ..col.column import Column +from ..constants.concrete_data import MAX_CONCRETE_STRAIN +from ..constants.rebar_data import STEEL_E, BarSize, rebar_area +from ..pmm_search.load_combo import LoadCombination + + +def getCalculatedColumnProps(bar_size: BarSize): + A_b = Calculation("Ab", rebar_area(bar_size)) + + E_s = Calculation("E_s", STEEL_E) + e_c = Calculation( + "\\epsilon_u", + MAX_CONCRETE_STRAIN, + ) + + return {"A_b": A_b, "E_s": E_s, "e_c": e_c} + + +@pytest.fixture +def example_col(): + bar_size: BarSize = "#5" + calc_props = getCalculatedColumnProps(bar_size) + + return Column( + Input("w", 24), + Input("h", 18), + Input("bar_size", bar_size), + Input("cover", 1.5), + Input("nx", 5), + Input("ny", 4), + Input("f'_c", 8000), + Input("f_y", 80), + False, + False, + calc_props["A_b"], + calc_props["E_s"], + calc_props["e_c"], + ) + + +@pytest.fixture +def loads(): + # for each load case: P, Mx, My, and whether the calc should be shown + loads = [[300, 100, 200, True], [-100, 50, -60, False], [11500, 300, -300, False]] + return [LoadCombination(i + 1, *load) for i, load in enumerate(loads)] + + +@pytest.fixture +def example_col2(): + bar_size: BarSize = "#8" + calc_props = getCalculatedColumnProps(bar_size) + + return Column( + Input("w", 16), + Input("h", 20), + Input("bar_size", bar_size), + Input("cover", 2.5), + Input("nx", 3), + Input("ny", 4), + Input("f'_c", 6000), + Input("f_y", 60), + True, + False, + calc_props["A_b"], + calc_props["E_s"], + calc_props["e_c"], + ) + + +@pytest.fixture +def example_col3(): + bar_size: BarSize = "#8" + calc_props = getCalculatedColumnProps(bar_size) + + return Column( + Input("w", 24), + Input("h", 36), + Input("bar_size", bar_size), + Input("cover", 1.5), + Input("nx", 3), + Input("ny", 4), + Input("f'_c", 4000), + Input("f_y", 40), + False, + False, + calc_props["A_b"], + calc_props["E_s"], + calc_props["e_c"], + ) diff --git a/examples/conc_col_pmm/tests/test_calculation.py b/examples/conc_col_pmm/tests/test_calculation.py new file mode 100644 index 0000000..87885bb --- /dev/null +++ b/examples/conc_col_pmm/tests/test_calculation.py @@ -0,0 +1,88 @@ +import matplotlib +import pytest + +from efficalc import Calculation, clear_saved_objects +from efficalc.calculation_runner import CalculationRunner + +from ..calc_document.calculation import calculation +from ..calc_document.column_inputs import ColumnInputs + +matplotlib.use("Agg") # Use a non-interactive backend + + +@pytest.fixture +def common_setup_teardown(): + # Set up a sample number + yield None # Provide the data to the test + # Teardown: Clean up resources (if any) after the test + clear_saved_objects() + + +def get_calc_by_name(all_items, name): + for item in all_items: + if isinstance(item, Calculation) and item.name == name: + return item + + +def assert_calc_value(calc: Calculation, expected: float): + assert calc.result() == pytest.approx(expected, abs=0.001) + + +def test_calc_with_defaults(common_setup_teardown): + runner = CalculationRunner(calculation) + all_obj = runner.calculate_all_items() + + ppn = get_calc_by_name(all_obj, "{\\phi}P_n") + pmx = get_calc_by_name(all_obj, "{\\phi}M_{nx}") + pmy = get_calc_by_name(all_obj, "{\\phi}M_{ny}") + dcr_mx = get_calc_by_name(all_obj, "DCR_{Mx}") + dcr_my = get_calc_by_name(all_obj, "DCR_{My}") + dcr_p = get_calc_by_name(all_obj, "DCR_{P}") + + assert_calc_value(ppn, 3579.613) + assert_calc_value(pmx, 238.932) + assert_calc_value(pmy, 119.412) + assert_calc_value(dcr_mx, 0.837060) + assert_calc_value(dcr_my, 0.837060) + assert_calc_value(dcr_p, 0.838079) + + +def test_calc_with_custom_load_case(common_setup_teardown): + loads = [[500, 400, 50, "yes"]] + runner = CalculationRunner(lambda: calculation(default_loads=loads)) + all_obj = runner.calculate_all_items() + + ppn = get_calc_by_name(all_obj, "{\\phi}P_n") + pmx = get_calc_by_name(all_obj, "{\\phi}M_{nx}") + pmy = get_calc_by_name(all_obj, "{\\phi}M_{ny}") + dcr_mx = get_calc_by_name(all_obj, "DCR_{Mx}") + dcr_my = get_calc_by_name(all_obj, "DCR_{My}") + dcr_p = get_calc_by_name(all_obj, "DCR_{P}") + + assert_calc_value(ppn, 2140.047) + assert_calc_value(pmx, 1712.806) + assert_calc_value(pmy, 214.404) + assert_calc_value(dcr_mx, 0.23353) + assert_calc_value(dcr_my, 0.23320) + assert_calc_value(dcr_p, 0.23364) + + +def test_calc_with_small_column(common_setup_teardown): + loads = [[18.22, 1.56, 3.03, "yes"]] + col = ColumnInputs(4, 6, "#4", 1, 2, 3, 4000, 40, True, True) + runner = CalculationRunner(lambda: calculation(default_loads=loads, col=col)) + all_obj = runner.calculate_all_items() + + ppn = get_calc_by_name(all_obj, "{\\phi}P_n") + pmx = get_calc_by_name(all_obj, "{\\phi}M_{nx}") + pmy = get_calc_by_name(all_obj, "{\\phi}M_{ny}") + dcr_mx = get_calc_by_name(all_obj, "DCR_{Mx}") + dcr_my = get_calc_by_name(all_obj, "DCR_{My}") + dcr_p = get_calc_by_name(all_obj, "DCR_{P}") + + assert_calc_value(ppn, 26.6350) + assert_calc_value(pmx, 2.28086) + assert_calc_value(pmy, 4.43000) + assert_calc_value(dcr_mx, 0.68395) + assert_calc_value(dcr_my, 0.68397) + assert_calc_value(dcr_p, 0.68406) diff --git a/examples/conc_col_pmm/tests/test_get_capacity.py b/examples/conc_col_pmm/tests/test_get_capacity.py new file mode 100644 index 0000000..45799d0 --- /dev/null +++ b/examples/conc_col_pmm/tests/test_get_capacity.py @@ -0,0 +1,14 @@ +from ..calc_document.plotting import get_capacity, pmm_mesh +from ..col.assign_max_min import calculate_axial_load_limits + + +def test_get_capacity(example_col, loads): + + col = example_col + axial_limits = calculate_axial_load_limits(col) + + # Retrieve the quarter PMM mesh, which has points + # in the format (Mx, My, P). + _, _, _, mesh = pmm_mesh.get_mesh(col, 48, 18, axial_limits) + + _ = get_capacity.get_capacity(mesh, loads[0]) diff --git a/examples/conc_col_pmm/tests/test_get_dcr_ecc.py b/examples/conc_col_pmm/tests/test_get_dcr_ecc.py new file mode 100644 index 0000000..106d39d --- /dev/null +++ b/examples/conc_col_pmm/tests/test_get_dcr_ecc.py @@ -0,0 +1,187 @@ +from ..col.assign_max_min import calculate_axial_load_limits +from ..pmm_search.ecc_search.get_dcr_ecc import get_dcr_ecc +from ..pmm_search.load_combo import LoadCombination + +""" +This test uses a set of load points and a given column as well as +reference values for the DCRs for those load points to check the +accuracy of the DCRs from this program. +""" +loads = [ + [-300, 50, 0], + [-300, 1000, 0], + [-200, 1000, 0], + [-100, 1000, 0], + [0, 1000, 0], + [100, 1000, 0], + [200, 1000, 0], + [300, 1000, 0], + [400, 1000, 0], + [500, 1000, 0], + [600, 1000, 0], + [700, 1000, 0], + [800, 1000, 0], + [900, 1000, 0], + [1000, 1000, 0], + [1100, 1000, 0], + [1200, 1000, 0], + [1300, 1000, 0], + [1400, 1000, 0], + [1500, 1000, 0], + [1600, 1000, 0], + [1700, 1000, 0], + [1800, 1000, 0], + [1900, 1000, 0], + [2000, 1000, 0], + [2100, 1000, 0], + [2200, 1000, 0], + [2200, 50, 0], + [-300, -50, 0], + [-300, -1000, 0], + [-200, -1000, 0], + [-100, -1000, 0], + [0, -1000, 0], + [100, -1000, 0], + [200, -1000, 0], + [300, -1000, 0], + [400, -1000, 0], + [500, -1000, 0], + [600, -1000, 0], + [700, -1000, 0], + [800, -1000, 0], + [900, -1000, 0], + [1000, -1000, 0], + [1100, -1000, 0], + [1200, -1000, 0], + [1300, -1000, 0], + [1400, -1000, 0], + [1500, -1000, 0], + [1600, -1000, 0], + [1700, -1000, 0], + [1800, -1000, 0], + [1900, -1000, 0], + [2000, -1000, 0], + [2100, -1000, 0], + [2200, -1000, 0], + [2200, -50, 0], + [-300, 0, 50], + [-300, 0, 1000], + [-200, 0, 1000], + [-100, 0, 1000], + [0, 0, 1000], + [100, 0, 1000], + [200, 0, 1000], + [300, 0, 1000], + [400, 0, 1000], + [500, 0, 1000], + [600, 0, 1000], + [700, 0, 1000], + [800, 0, 1000], + [900, 0, 1000], + [1000, 0, 1000], + [1100, 0, 1000], + [1200, 0, 1000], + [1300, 0, 1000], + [1400, 0, 1000], + [1500, 0, 1000], + [1600, 0, 1000], + [2000, 0, 0], +] + +dcrs = [ + 1.18111455, + 3.5800972, + 3.24327159, + 2.90662146, + 2.56997156, + 2.23332143, + 1.89667141, + 1.60175514, + 1.34322953, + 1.148969, + 1.01194346, + 0.911275268, + 0.9276003, + 0.9705417, + 1.01678479, + 1.06635, + 1.11907816, + 1.15417635, + 1.184222, + 1.217052, + 1.25151849, + 1.28756726, + 1.324892, + 1.36372554, + 1.40255916, + 1.44369924, + 1.48486662, + 1.31116092, + 1.18111455, + 3.5800972, + 3.24327159, + 2.90662146, + 2.56997156, + 2.23332119, + 1.896671, + 1.601755, + 1.34322941, + 1.148969, + 1.01194358, + 0.9112752, + 0.9276003, + 0.9705417, + 1.01678479, + 1.06635, + 1.11907816, + 1.15417635, + 1.184222, + 1.217052, + 1.25151849, + 1.28756726, + 1.324892, + 1.36372554, + 1.40255916, + 1.44369924, + 1.48486662, + 1.31116092, + 1.23884523, + 4.9428153, + 4.608767, + 4.274719, + 3.94067168, + 3.60662413, + 3.27257633, + 2.93852878, + 2.62903023, + 2.32137275, + 2.022534, + 1.78518534, + 1.59556031, + 1.4507966, + 1.35867274, + 1.35752714, + 1.39450872, + 1.44129765, + 1.49240327, + 1.54428363, + 1.59789062, + 1.19196439, +] + +# the error tolerance for this test +tol = 1e-2 + + +def test_get_dcr(example_col3): + col = example_col3 + axial_limits = calculate_axial_load_limits(col) + + for i in range(78): + load = LoadCombination(i + 1, *(loads[i]), False) + dcr = get_dcr_ecc(col, load, axial_limits) + + if dcr > 0: + assert abs((dcr - dcrs[i]) / dcr) < tol + else: + assert dcrs[0] < tol diff --git a/examples/conc_col_pmm/tests/test_pmm_plotter_plotly.py b/examples/conc_col_pmm/tests/test_pmm_plotter_plotly.py new file mode 100644 index 0000000..8649e2c --- /dev/null +++ b/examples/conc_col_pmm/tests/test_pmm_plotter_plotly.py @@ -0,0 +1,20 @@ +from ..calc_document.plotting import get_pmm_data, pmm_plotter_plotly +from ..col import assign_max_min +from ..pmm_search.load_combo import LoadCombination + + +# This test checks for runtime errors +def test_pmm_plotter_plotly(example_col, loads): + axial_limits = assign_max_min.calculate_axial_load_limits(example_col) + + # for each load case: P, Mx, My, and whether the calc should be shown + load_data = [ + [300, 100, 200, True], + [-100, 50, -60, False], + [1500, 300, -300, False], + ] + loads = [LoadCombination(i, *load) for i, load in enumerate(load_data)] + + pmm_data = get_pmm_data.get_pmm_data(example_col, 36, 12, loads, axial_limits) + + _ = pmm_plotter_plotly.plot(pmm_data) diff --git a/examples/conc_col_pmm/tests/test_point_plotter.py b/examples/conc_col_pmm/tests/test_point_plotter.py new file mode 100644 index 0000000..78f1dc1 --- /dev/null +++ b/examples/conc_col_pmm/tests/test_point_plotter.py @@ -0,0 +1,13 @@ +from ..calc_document.plotting import get_capacity, pmm_mesh, point_plotter +from ..col.assign_max_min import calculate_axial_load_limits + +# This test checks for runtime errors + + +def test_point_plotter(example_col3, loads): + + axial_limits = calculate_axial_load_limits(example_col3) + _, _, _, mesh = pmm_mesh.get_mesh(example_col3, 36, 12, axial_limits) + + capacity = get_capacity.get_capacity(mesh, loads[0]) + _ = point_plotter.plot(capacity, loads[0], False) diff --git a/examples/conc_col_pmm/tests/test_point_search_ecc.py b/examples/conc_col_pmm/tests/test_point_search_ecc.py new file mode 100644 index 0000000..5197a82 --- /dev/null +++ b/examples/conc_col_pmm/tests/test_point_search_ecc.py @@ -0,0 +1,34 @@ +import math + +from ..col.assign_max_min import calculate_axial_load_limits +from ..pmm_search.ecc_search.point_search_ecc import search + +search_tol = 1.5e-3 # the tolerance for error in the points found +ceil_tol = 1e-2 # how close to limits to go + +""" +The function below chooses an arbitrary initial guess and then +for a range of target points (where the target lambda and load +are both allowed to vary) checks whether the eccentricity search +algorithm converges. +""" + + +def test_search(example_col): + axial_limits = calculate_axial_load_limits(example_col) + guess = [-0.7, 30] + lambda_change = math.pi / 12 + ecc_change = math.pi / 12 + lambda_target = 0 + while lambda_target < math.pi / 2 - ceil_tol: + ecc_target = -math.pi / 2 + ecc_change + while ecc_target < math.pi / 2 - ceil_tol: + target = [lambda_target, ecc_target] + Mx, My, P, guess = search(example_col, target, guess, axial_limits) + lambda_found = math.atan2(My, Mx) + Mxy_found = math.sqrt(Mx**2 + My**2) + ecc_found = math.atan2(P, Mxy_found) + assert abs(lambda_found - lambda_target) < search_tol + assert abs(ecc_found - ecc_target) < search_tol + ecc_target += ecc_change + lambda_target += lambda_change diff --git a/examples/conc_col_pmm/tests/test_point_search_load.py b/examples/conc_col_pmm/tests/test_point_search_load.py new file mode 100644 index 0000000..5e4048d --- /dev/null +++ b/examples/conc_col_pmm/tests/test_point_search_load.py @@ -0,0 +1,32 @@ +import math + +from ..col.assign_max_min import calculate_axial_load_limits +from ..pmm_search.load_search.point_search_load import search + +search_tol = 1.5e-3 # the tolerance for error in the points found +ceil_tol = 1e-2 # how close to limits to go + +""" +The function below chooses an arbitrary initial guess and then +for a range of target points (where the target lambda and load +are both allowed to vary) checks whether the load search +algorithm converges. +""" + + +def test_search(example_col): + axial_limits = calculate_axial_load_limits(example_col) + guess = [-0.7, 30] + lambda_change = math.pi / 12 + load_change = axial_limits.load_span / 10 + lambda_target = 0 + while lambda_target < math.pi / 2 - ceil_tol: + load_target = axial_limits.min_phi_pn + load_change + while load_target < axial_limits.max_phi_pn - ceil_tol: + target = [lambda_target, load_target] + Mx, My, P, _, _ = search(example_col, target, guess, axial_limits) + lambda_found = math.atan2(My, Mx) + assert abs(lambda_found - lambda_target) < search_tol + assert (abs(P - load_target) / axial_limits.load_span) < search_tol + load_target += load_change + lambda_target += lambda_change diff --git a/examples/conc_col_pmm/tests/test_triangles.py b/examples/conc_col_pmm/tests/test_triangles.py new file mode 100644 index 0000000..048e89f --- /dev/null +++ b/examples/conc_col_pmm/tests/test_triangles.py @@ -0,0 +1,11 @@ +from ..struct_analysis.triangles import * + +tri1 = ((0, 0), (3, 0), (0, 6)) + + +def test_triangle_area(): + assert triangle_area(tri1[0], tri1[1], tri1[2]) == 9 + + +def test_triangle_centroid(): + assert triangle_centroid(tri1[0], tri1[1], tri1[2]) == (1, 2) diff --git a/examples/conc_col_pmm/tests/test_try_axis.py b/examples/conc_col_pmm/tests/test_try_axis.py new file mode 100644 index 0000000..11d8f84 --- /dev/null +++ b/examples/conc_col_pmm/tests/test_try_axis.py @@ -0,0 +1,33 @@ +import math + +from ..col.assign_max_min import calculate_axial_load_limits +from ..struct_analysis.try_axis import try_axis + +# this test is based on the reference calculation by SP Column which can be found at the following link: +# https://structurepoint.org/publication/pdf/Biaxial-Bending-Interaction-Diagrams-for-Rectangular-Reinforced-Concrete-Column-Design-ACI-318-19.pdf +# page 11 of the page above indicates the values of theta and c being tried as well as the values +# for factored Mx, My, and P using the exact capacity method. + +# The reference P, Mx, My are below, these are factored +reference_forces = (283.72, 214.29, 133.83) + +# the error tolerance for this test +tol = 1e-3 + + +def test_try_axis(example_col2): + axial_limits = calculate_axial_load_limits(example_col2) + + (lambda1, ecc, pn, phi_pn_not_limited, phi_mn_xy) = try_axis( + example_col2, -43.9 * math.pi / 180, 12.5, axial_limits + ) + phi_pn = min(axial_limits.max_phi_pn, phi_pn_not_limited) + + found = ( + phi_pn, + phi_mn_xy * math.cos(lambda1), + phi_mn_xy * math.sin(lambda1), + ) + for i in range(3): + ref = reference_forces[i] + assert abs((found[i] - ref) / ref) < tol diff --git a/examples/conc_col_pmm/tests/test_try_axis_document.py b/examples/conc_col_pmm/tests/test_try_axis_document.py new file mode 100644 index 0000000..82b3866 --- /dev/null +++ b/examples/conc_col_pmm/tests/test_try_axis_document.py @@ -0,0 +1,31 @@ +import math + +from ..calc_document.try_axis_document import try_axis_document +from ..col.assign_max_min import calculate_axial_load_limits + +# this test is based on the reference calculation by SP Column which can be found at the following link: +# https://structurepoint.org/publication/pdf/Biaxial-Bending-Interaction-Diagrams-for-Rectangular-Reinforced-Concrete-Column-Design-ACI-318-19.pdf +# page 11 of the page above indicates the values of theta and c being tried as well as the values +# for factored Mx, My, and P using the exact capacity method. + +# The reference P, Mx, My are below, these are factored +reference_forces = (283.72, 214.29, 133.83) + +# the error tolerance for this test +tol = 1e-3 + + +def test_try_axis(example_col2): + axial_limits = calculate_axial_load_limits(example_col2) + capacities = try_axis_document( + example_col2, axial_limits, -43.9 * math.pi / 180, 12.5 + ) + found = ( + capacities.P.result(), + capacities.Mx.result(), + capacities.My.result(), + ) + print("found", found) + for i in range(3): + ref = reference_forces[i] + assert abs((found[i] - ref) / ref) < tol diff --git a/examples/conc_col_pmm/tests/visual_tests/__init__.py b/examples/conc_col_pmm/tests/visual_tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/examples/conc_col_pmm/tests/visual_tests/visual_test_document_wrapper.py b/examples/conc_col_pmm/tests/visual_tests/visual_test_document_wrapper.py new file mode 100644 index 0000000..09c0bd0 --- /dev/null +++ b/examples/conc_col_pmm/tests/visual_tests/visual_test_document_wrapper.py @@ -0,0 +1,29 @@ +import os +import sys + +# Now you can import from conc_col_pmm +from ...calc_document.document_wrapper import run + +# "w", "h", "bar_size", "bar_cover", "bars_x", "bars_y", "fc", "fy", "cover_type", "transverse_type", +col_data = [24, 18, "#6", 2, 5, 2, 8000, 60, "Edge", "Spiral"] + +# for each load case: P, Mx, My, and whether the calc should be shown +# Note that these load cases currently do not override the defaults +loads = [ + [300, 100, 200, "yes"], + [-100, 50, -60, "no"], + [11500, 300, -300, "no"], + [0, 200, 0, "yes"], + [0, 0, 200, "yes"], +] + +# calc_report_example1 +# col_data = [24, 18, "#6", 1.5, 5, 4, 8000, 60, "Edge", "Tied"] +# loads = [[1400, -300, 100, True]] + +# calc_report_example2 +# col_data = [24, 36, "#8", 2, 6, 8, 8000, 60, "Edge", "Tied"] +# loads = [[3000, -200, 100, True]] + +if __name__ == "__main__": + run(True, col_data, loads) diff --git a/examples/conc_col_pmm/tests/visual_tests/visual_test_pmm_plotter_plotly.py b/examples/conc_col_pmm/tests/visual_tests/visual_test_pmm_plotter_plotly.py new file mode 100644 index 0000000..f3ae628 --- /dev/null +++ b/examples/conc_col_pmm/tests/visual_tests/visual_test_pmm_plotter_plotly.py @@ -0,0 +1,50 @@ +from efficalc import Input + +from ...calc_document.plotting import get_pmm_data, pmm_plotter_plotly +from ...col import assign_max_min +from ...col.column import Column +from ...constants.rebar_data import BarSize +from ...pmm_search.load_combo import LoadCombination +from ...tests.conftest import getCalculatedColumnProps + +# TODO: make this use the main calc callsite and get the plotly data from there + + +def example_col(): + bar_size: BarSize = "#6" + calc_props = getCalculatedColumnProps(bar_size) + + return Column( + Input("w", 24), + Input("h", 18), + Input("bar_size", bar_size), + Input("cover", 2), + Input("nx", 5), + Input("ny", 2), + Input("f'_c", 8000), + Input("f_y", 60), + False, + True, + calc_props["A_b"], + calc_props["E_s"], + calc_props["e_c"], + ) + + +if __name__ == "__main__": + col = example_col() + axial_limits = assign_max_min.calculate_axial_load_limits(col) + + # for each load case: P, Mx, My, and whether the calc should be shown + load_data = [ + [300, 100, 200, True], + [-100, 50, -60, False], + [1500, 300, -300, False], + ] + loads = [LoadCombination(i, *load) for i, load in enumerate(load_data)] + + pmm_data = get_pmm_data.get_pmm_data(col, 36, 12, loads, axial_limits) + + fig = pmm_plotter_plotly.plot(pmm_data) + + fig.show() diff --git a/examples/conc_col_pmm/tests/visual_tests/visual_test_point_plotter.py b/examples/conc_col_pmm/tests/visual_tests/visual_test_point_plotter.py new file mode 100644 index 0000000..4704d0a --- /dev/null +++ b/examples/conc_col_pmm/tests/visual_tests/visual_test_point_plotter.py @@ -0,0 +1,53 @@ +import matplotlib.pyplot as plt + +from efficalc import Input + +from ...calc_document.plotting import get_capacity, pmm_mesh, point_plotter +from ...col import assign_max_min +from ...col.column import Column +from ...constants.rebar_data import BarSize +from ...pmm_search.load_combo import LoadCombination +from ...tests.conftest import getCalculatedColumnProps + + +def example_col(): + bar_size: BarSize = "#5" + calc_props = getCalculatedColumnProps(bar_size) + + return Column( + Input("w", 24), + Input("h", 18), + Input("bar_size", bar_size), + Input("cover", 1.5), + Input("nx", 5), + Input("ny", 4), + Input("f'_c", 8000), + Input("f_y", 80), + False, + False, + calc_props["A_b"], + calc_props["E_s"], + calc_props["e_c"], + ) + + +if __name__ == "__main__": + col = example_col() + axial_limits = assign_max_min.calculate_axial_load_limits(col) + + # for each load case: P, Mx, My, and whether the calc should be shown + load_data = [ + [300, 100, 200, True], + [-100, 50, -60, False], + [11500, 300, -300, False], + ] + loads = [LoadCombination(i, *load) for i, load in enumerate(load_data)] + + _, _, _, mesh = pmm_mesh.get_mesh(col, 36, 12, axial_limits) + + capacity = get_capacity.get_capacity(mesh, loads[0]) + fig = point_plotter.plot(capacity, loads[0], False) + + plt.show() + + # plt.savefig("test_plot.png") diff --git a/requirements.txt b/requirements.txt index ce76d93..1a522ac 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ latexexpr_efficalc==0.5.3 pylatexenc==2.10 -pytest==8.0.2 +pytest==8.0.2 \ No newline at end of file