Experiment/results-object

[talledodiego]

This is just an experiment of what I was trying before our discussion @mortenengen.

Outdated, new developments later

The current implementation is really easy even if works only with fibers it does not requires additional methods to geometry since it deals everything from triangulation data already available with fibers.

But this basically permits doing something simple like:

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
y, z, A, eps = res.get_fibers_strains()
y, z, A, sig = res.get_fibers_stresses()
axs[0].scatter(y, z, c=eps, s=10)
axs[0].axis('equal')
axs[1].scatter(y, z, c=sig, s=10)
axs[1].axis('equal')

Obtaining:

Or even more elaborated things like:

fig, axs = plt.subplots(1, 2, figsize=(12, 5))
y, z, A, eps = res.get_fibers_strains()
yc, zc, Ac, sigc = res.get_fibers_stresses(0)
ys, zs, As, sigs = res.get_fibers_stresses(1)
yp, zp, Ap, sigp = res.get_fibers_stresses(2)
sc = axs[0].scatter(y, z, c=eps, s=10)
plt.colorbar(sc, ax=axs[0], label='Strain')

axs[0].axis('equal')
sc = axs[1].scatter(yc, zc, c=sigc, s=10)
plt.colorbar(sc, ax=axs[1], label=r'$\sigma_c$ [MPa]')
reinf_min = min(min(sigs), min(sigp))
reinf_max = max(max(sigs), max(sigp))
sc = axs[1].scatter(yp, zp, c=sigp, s=80, edgecolors='k', cmap='cool', vmin=reinf_min, vmax=reinf_max)
axs[1].scatter(ys, zs, c=sigs, s=30, edgecolors='k', cmap=sc.cmap)
plt.colorbar(sc, ax=axs[1], label=r'$\sigma_s$ [MPa]')

axs[1].axis('equal')

plt.tight_layout()

Obtaining as an example this:

The relation between coordinates and strains at each coordinate is:

$$\varepsilon_a - k_z y + k_y z$$

That in matrix form is:

$$ \varepsilon = \begin{bmatrix} 1 & -y & z \end{bmatrix} \begin{bmatrix} \varepsilon_a \ k_z \ k_y \end{bmatrix}^T $$

I find somehow strange (but I can think more about it) to have this vector

$$ \begin{bmatrix} 1 & -y & z \end{bmatrix} $$

as a return of the geometry since a "geometry" is nothing till when it is not inside the section, and it is when it is in the section that has that specific reference system and therefore that the vector

$$ \begin{bmatrix} \varepsilon_a \ k_z \ k_y \end{bmatrix}^T $$

assumes that meaning...

I don't know if it is just me, but I find more logic having this stuff in the level of sections (or in this case results), not geometries.

New implementations

I created the following new feature in the geometries that permits us to better code the methods seen above and to make them work not only for fiber integrator, but also for marin integrator.
Moreover having such dense points within the polygon it is possible to add a probe method passing coordinates of a point (y,z) and returning the strain and/or stress responde of the nearest point to the input one.

Get random points within polygon

Meanwhile I did a first fast implementation to get a bunch of random points uniformly distributed from the SurfaceGeometry.

I considered two possible algorithms, but I choose the fastest one.

For instance for this polygon:

we get (these are 1000 points):

Or for this geometry:

poly = RectangularGeometry(width=600, height=400, material=concrete)
hole = RectangularGeometry(width=200, height=100, material=concrete)
poly -= hole
poly

with this simple command:

x,y = poly.random_points_within(num_points=10000)
plt.plot(x,y,'o')
plt.axis('equal');

we obtain instantly:

This is the path for permitting not using triangulation data (like in my current test I had) making it possible to work with both fiber or marin integrators.

Get point coordinates as a numpy array from compound geometry

I also added a method for returning a numpy array from a CompoundGeometry. The method by default returns all point geometries. If the optional group_label parameter is passed, this is used as a filter:

# Return all point geometries coordinates
geo.get_coordinates_point_geometries()

# Return all point geometries coordinates of group "reinforcement"
geo.get_coordinates_point_geometries(group_label="prestress_reinforcement")

If the filter is not found the array will be empty.

New detailed results

As an attempt I created a SectionDetailedResultState that basically stores strains and stress for all the geometry for a given state (i.e. a set of deformations / forces). For instance when they come from a computation of ultimate strength it will be the failure state, but during a moment-curvature each state corresponds to one point of the curve.

For now as an attempt for testing there is a method for ceating this datastructure that remains stored (up to when we don't want to change state). In this way we can inspect (or plot) for each part of the geometry strains and / or stresses.
By default the point geometries are grouped by group_label.

This example shows the idea.

from structuralcodes import set_design_code
from structuralcodes.materials.concrete import create_concrete
from structuralcodes.materials.reinforcement import create_reinforcement
from structuralcodes.geometry import RectangularGeometry, add_reinforcement_line, add_reinforcement, SurfaceGeometry
from structuralcodes.sections import GenericSection

from shapely.geometry import Polygon
import matplotlib.pyplot as plt
import numpy as np

set_design_code('ec2_2004')

concrete = create_concrete(45, alpha_cc=0.85)
print(f'Concrete fcd = {concrete.fcd():.2f} MPa')

reinforcement = create_reinforcement(450, 200000, 450, 0.075)
print(f'Reinforcement fyd = {reinforcement.fyd():.2f} MPa')

# Create the prestressing reinforcement with initial strain
prestress_reinf = create_reinforcement(fyk=1620, Es=200000, ftk=1620, epsuk=0.02, initial_strain=0.0054)
print(f'Initial stress = {prestress_reinf.initial_stress:.2f} MPa') 
print(f'Initial strain = {prestress_reinf.initial_strain*1e3:.2f} ‰') 

poly = Polygon([
    (250, 0), 
    (250, 140), 
    (70, 200), 
    (70, 920), 
    (250, 980), 
    (250, 1120), 
    (-250, 1120), 
    (-250, 980), 
    (-70, 920), 
    (-70, 200), 
    (-250, 140), 
    (-250, 0),
])

geo = SurfaceGeometry(poly=poly, material=concrete)
for y, n in zip([40, 110, 1010, 1080], [4, 2, 2, 4]):
    geo = add_reinforcement_line(
        geo=geo,
        coords_i=(-210, y), 
        coords_j=(210, y),
        diameter=10,
        n=n,
        material=reinforcement,
        group_label='reinforcement'
    )
area = 840
d = (area * 4 / np.pi)**0.5
geo = add_reinforcement_line(
    geo=geo,
    coords_i=(0, 70),
    coords_j=(0, 210),
    diameter=d,
    n=2,
    material=prestress_reinf,
    group_label='prestress_reinforcement'
)
geo

section = GenericSection(geometry=geo, integrator='fiber', mesh_size=0.001)
res = section.section_calculator.calculate_bending_strength()
print(f'Bending strength: {res.m_y*1e-6:.2f} kNm')


res.create_detailed_result()

now res contains a detailed_result object which contains surface_data and point_data. The plotting is slightly more complicated than what I have done above for fibers, but it permits way much liberty:

# Create a plot of strains
y = np.empty(0)
z = np.empty(0)
eps = np.empty(0)
for sd in res.detailed_result.surface_data:
    y = np.append(y, sd['y'])
    z = np.append(z, sd['z'])
    eps = np.append(eps, sd['strain'])
sc = plt.scatter(y, z, c=eps)
plt.colorbar(sc, label='Strain')
plt.axis('equal');

or

# Create a plot of strains
y = np.empty(0)
z = np.empty(0)
sig = np.empty(0)
for sd in res.detailed_result.surface_data:
    y = np.append(y, sd['y'])
    z = np.append(z, sd['z'])
    sig = np.append(sig, sd['stress'])
sc = plt.scatter(y, z, c=sig)
plt.colorbar(sc, label='Stress [MPa]')
plt.axis('equal');

obtaining for example:

But one could do a different color map for each surface if we have more SurfaceGeometries.

In the same way one could deal with point_data, plotting (or asking for) only certain group_label or all. For instance adding prestress reinforcement and regular reinforcement (by group_label) with different color maps and bars:

# Create a plot of strains
y = np.empty(0)
z = np.empty(0)
sig = np.empty(0)
for sd in res.detailed_result.surface_data:
    y = np.append(y, sd['y'])
    z = np.append(z, sd['z'])
    sig = np.append(sig, sd['stress'])
sc = plt.scatter(y, z, c=sig)
# show reinforcemnet results
coords = res.detailed_result.point_data['reinforcement']['coords']
sig = res.detailed_result.point_data['reinforcement']['stress']
sc_r = plt.scatter(coords[:, 0], coords[:, 1], c=sig, cmap = 'cool', edgecolors='k', s=20)
#show prestress reinforcment results
coords = res.detailed_result.point_data['prestress_reinforcement']['coords']
sig = res.detailed_result.point_data['prestress_reinforcement']['stress']
sc_p = plt.scatter(coords[:, 0], coords[:, 1], c=sig, cmap = 'GnBu', edgecolors='k', s=60)
plt.colorbar(sc, label=r'$\sigma_c$ [MPa]')
plt.colorbar(sc_r, label=r'$\sigma_{r}$ [MPa]')
plt.colorbar(sc_p, label=r'$\sigma_{p}$ [MPa]')

plt.axis('equal')
plt.tight_layout();

obtaining:

Finally this detailed_result object could contain a probe method passing the coordinates and returning strains and stresses. Some filtering utilities could be included in the probe for instance for inspecting only surfaces, surface with a name or label, points (also with names or labels).

New detailed results for moment-curvature

I also experimented with moment curvature, permitting to advance from DetailState to DetailState. In this way one can advance from one step to next with res.next_step() or go back with res.previous_step(); it is also possible to go to a specific step with res.set_step(). Each one of these commands prepare the detailed results objects that therefore can be used for plotting, querying etc.
Playing with such functions for instance one can create an animation of strain or stress distribution during the moment-curvature analysis, like the following animation, where it is visible the suddn shift in neautral axis after yielding of bars.

[mortenengen]

Could this portion of the code be refactored out to avoid the code repetition below? We could e.g. create an internal method _create_detailed_result which sets self.detailed_result based on self.current_step, so for each of the methods create_detailed_result, next_step, previous_step and set_step we can simply set self.current_step and then call the internal method.

[mortenengen]

It would be great if we could store the surface data in a hierarchy similar to the point data, e.g. as a dict where the key is the group_label and the value is the dict with data.

Alternatively, and consistently for surface and point data, we could store the data in a manner that is easily input to a dataframe. For this structure, I guess one option could be to store a dict with the following keys: x, y, strain, stress, group_label, and each value is an array or a list of values for each individual point. So one dict for surface data and one for point data. This would open for really smooth filtering and plotting downstream.

[mortenengen]

There are still lines that are not covered in testing. Could we perhaps add a test where we for example create detailed results for a moment curvature relation, and for each set of detailed results, integrate over the results to check if the detailed results add up to the moment, axial force, curvature and axial strain of the moment curvature relation?

[mortenengen]

Thanks for this valuable feature addition! I have left a couple of comments for you to consider during finalizing 😃