Issues with creating geometry from points.

[EngDM]

Hey all,

In the documentation under section 3.2.1 it mentions that if you don't have a complex shape that you do not need to provide facets, however when running the code it seems that I do. Am I missing something? Also, what constitutes a facet in the input? I tried deciphering the correlation from the values provided in the example to the shape drawn, but I can't see how the coordinates provided indicate a face of the shape.

Below is my code:

from sectionproperties.pre.geometry import Geometry, CompoundGeometry
from sectionproperties.analysis.section import Section
import math as m

#INPUT YOUR VARIABLES (2 decimal points)
a = float(input('Leg Length a: '))
t = float(input('Member Thickness t: '))
r = float(input('Angle α: '))
c = float(input('Bottom Length c: '))
d = float(input('Joist Depth d: '))
l = float(input('Joist Length l: '))

To replicate my shape input the following:

x = m.sin(m.radians((90-r)/2))*(t/(m.sin(m.radians(90-(90-r)/2))))
points = [
    [0,0],
    [a+x,0],
    [(l-c)/2+x,t-d],
    [(l+c)/2-x,t-d],
    [l-a-x,0],
    [l,0],
    [l,-t],
    [l-a,-t],
    [(l+c)/2,-d],
    [(l-c)/2,-d],
    [a,-t],
    [0,-t]
holes = []
facets = []
control_points = [[0,0]]
]

EDIT: I have got it to work by inputting completely ambiguous values for facets similar to the example, but I haven't the slightest clue why it is working.
EDIT2: I misinterpreted the documentation in that a facet was a coordinate, but instead it is a list of the points as if we were connecting the dots. My assumption is that the integers represent the tuples inputted into the points list in the order that they appear.

facets = [
[0, 1],
[1, 2],
[2, 3],
[3, 4],
[4, 5],
[5, 0],
[6, 7],
[7, 8],
[8, 9],
[9, 6],
[10,5],
[11,6],
]

facets = [
    [0, 1],
    [1, 2],
    [2, 3],
    [3, 4],
    [4, 5],
    [5, 6],
    [6, 7],
    [7, 8],
    [8, 9],
    [9, 10],
    [10,11],
    [11,12],
]

geometry = Geometry.from_points(points, facets, control_points)
geometry.create_mesh([1])

I appreciate the input and for your patience with a Python novice.

[EngDM]

I thought I'd keep this in a separate post as it is more of an inquiry.
I'm trying to align the reinforcing shown (pink circle) with some typical shapes, however the align_to bottom puts the circle way far out.

#Reinforcing Shapes
pointsr1 = circular_section(d=12.7, n=17).align_to(geometry, on="bottom").align_to(geometry, on="left")
pointsr2 = circular_section(d=12.7, n=17).align_to(geometry, on="bottom").align_to(geometry, on="right")

#Top Chord Shape
geometry = geometry + pointsr1

Update:

I've developed an equation relating the corner point and the point in which the circle can be centered on so that the corners are perfectly aligned. I'll post the code tomorrow, but the equation is as follows:

$d=r*(cos(α)-(1-sin(α))*tan(α))$

Where d is the horizontal distance you need to move, and α is the angle between side of the vertical leg and the Y axis vertical plane and r is the radius of your circle. Whether there is a way to simplify this equation, I'm not sure. But it looks clean.

[Agent6-6-6]

The align_to method only aligns to the bounding box of your geometry as noted in the documentation. So what you're seeing is the intended behaviour. I don't think there is a method built in to align the edges of geometries that you're looking for without going back to utilising some of the underlying shapely methods to determine the distance you need to move.

Alternatively if this is something you're only doing once and you know the coordinates of the underside of that shape then you can just move it manually.

[EngDM]

Ah okay, I'll play around with perhaps aligning to an x,y. Perhaps I can determine a relationship between the corner node where I want the circle and where the the edges align.

[robbievanleeuwen]

First question - the facets do represent the indices of the points, as described in the docs "A list of (start, end) indexes of vertices defining the edges of the section geometry." Unfortunately it seems like the documented optional facet list was never implemented, I'll raise an issue and get this fixed. In the meantime a workaround is the following:

from shapely.geometry import Polygon
from sectionproperties.pre.geometry import Geometry

points = [...]  # list of points

geometry = Geometry(Polygon(points))

[EngDM]

I also found in there that if you use the create_facets class on the list of points it does it for you, so it's like you are skipping the facets anyways.

[EngDM]

The final code is below, however for whatever reason the shear center calculation is inconsistent and sometimes it shoots it way off 200 units away, and other times it is where I expect it to be. This happens even just re-running identical code. There should be no overlapping sections so I'm not sure why this would happen.

Some runs it gives an error about point discontinuity (without the reinforcing added), and other times it does not.

Anyways, here it is for a preliminary, I plan to get it to run some calculations based on CSA S16 to calculate compressive resistance and moment resistance eventually, and loop thru different round bar sizes. But that beast is for another day.

from sectionproperties.pre.geometry import Geometry, CompoundGeometry
from sectionproperties.analysis.section import Section
from sectionproperties.pre.library.primitive_sections import circular_section
from IPython.display import Image
import sectionproperties.pre.geometry
import easygui
import math as m

a = float(easygui.enterbox('Leg Length a: ')) 
t = float(easygui.enterbox('Member Thickness t: '))
c = float(easygui.enterbox('Bottom Length c: '))
d = float(easygui.enterbox('Joist Depth d: '))
l = float(easygui.enterbox('Joist Length l: '))
r = abs(m.degrees(m.atan(((l-c)/2-a)/(t-d))))

x = m.sin(m.radians((90-r)/2))*(t/(m.sin(m.radians(90-(90-r)/2))))
points = [
    [0,0],
    [a+x,0],
    [(l-c)/2+x,t-d],
    [(l+c)/2-x,t-d],
    [l-a-x,0],
    [l,0],
    [l,-t],
    [l-a,-t],
    [(l+c)/2,-d],
    [(l-c)/2,-d],
    [a,-t],
    [0,-t],
]
facets = sectionproperties.pre.geometry.create_facets(points, True)
holes = []
control_points = [[0,0]]

#Top Chord Shape
geometry = Geometry.from_points(points, facets, control_points)

#Diameter of reinforcing (0 if none)
diam = float(easygui.enterbox('Enter Reinforcing Diameter: '))
#Offset of round bar from corner
y_off = diam/2
angle = abs(m.atan(((l-c)/2-a)/(t-d)))
x_off = y_off*(m.cos(angle)-(1-m.sin(angle))*m.tan(angle))

#Reinforcing Shapes
if diam>0:
    pointsr1 = circular_section(d=diam, n=20)
    pointsr1 = pointsr1.shift_section(x_offset = a-x_off, y_offset = -t-y_off)
    pointsr2 = circular_section(d=diam, n=20)
    pointsr2 = pointsr2.shift_section(x_offset = l-a+x_off, y_offset = -t-y_off)
else:
    pointsr1 = None
    pointsr2 = None
#Sum the geometry
geometry = geometry + pointsr1 + pointsr2
geometry = geometry.align_center()

#Plot and calculate properties
geometry.plot_geometry()
geometry.create_mesh([1])
Props = Section(geometry)
Props.calculate_geometric_properties()
Props.calculate_warping_properties()
Props.plot_mesh(materials=False)
Props.plot_centroids()
Props.display_results(fmt='.4f')

[Agent6-6-6]

@EngDM you can do that marginal overlap and trimming natively in section properties using techniques as per docs under advanced geometry creation.
https://sectionproperties.readthedocs.io/en/latest/rst/advanced_geom.html

It's a matter of ensuring the bars and the shape overlap enough (to not cause any issues as outlined in @robbievanleeuwen post which might require some experimentation), then create the geometry using:-

composite_shape = (bars - shape) + shape

This removes part of the bar that overlaps the other shape (like your manual trim), then just add the shape back to create the compound geometry that consists of the bars and the other shape.

[EngDM]

@Agent6-6-6 I like this approach, instead of moving the bars over by the exact geometrical amount, I'll set it to like 99% or whatever works to ensure overlap, maxing out at maybe 1mm or something TBD.

Edit:

I've implemented the above code for the composite shape, and for whatever reason now it does not want to run any mesh creation or geometry plotting, it runs the code indefinitely and I have to force stop it. Below is the shape with the nodes in common being shown.

This is the code that I ran on this (amended from the original code posted above)

if diam>0:
    pointsr1 = circular_section(d=diam, n=20)
    pointsr1 = pointsr1.shift_section(x_offset = a-x_off, y_offset = -t-y_off)
    
    
    pointsr2 = circular_section(d=diam, n=20)
    pointsr2 = pointsr2.shift_section(x_offset = l-a+x_off, y_offset = -t-y_off)
    
    nodes_in_common = (geometry - (pointsr1 + pointsr2)) | geometry
    
    analysis_geom = ((pointsr1 + pointsr2) - geometry) + nodes_in_common
   
    nodes_in_common.plot_geometry()
    
else:
    pointsr1 = None
    pointsr2 = None    
    analysis_geom = geometry.align_center()


analysis_geom.plot_geometry()
analysis_geom.create_mesh([10])
analysis_sec = cross_section.Section(analysis_geom)
analysis_sec.plot_mesh()
analysis_sec.calculate_geometric_properties()
analysis_sec.calculate_warping_properties()
analysis_sec.plot_centroids()
analysis_sec.display_results(fmt='.4f')

If I comment out everything below analysis_geom.plot_geometry() then it generates the combined shape, it just doesn't want to mesh it.

[EngDM]

@robbievanleeuwen Any idea why I would be getting the error "C:\Users~~~\AppData\Roaming\jupyterlab-desktop\jlab_server\lib\site-packages\sectionproperties\analysis\section.py:360: UserWarning:
The section geometry contains disjoint regions which is invalid for warping analysis.
Please revise your geometry to ensure there is connectivity between all regions.
Please see https://sectionproperties.readthedocs.io/en/latest/rst/analysis.html#warping-analysis for more information.
warnings.warn(" with the code above (without applying any reinforcing)? Does it think the reinforcing is a region even tho I input zero for diameter? This error is not consistent, only happens sometimes but can happen or not happen even with the same inputs.

[normanrichardson]

@EngDM
This is the code responsible for that warning:

section-properties/sectionproperties/analysis/section.py

Lines 355 to 365 in ca86307

	# check if any geometry are disjoint
	if isinstance(self.geometry, section_geometry.CompoundGeometry):
	polygons = [sec_geom.geom for sec_geom in self.geometry.geoms]
	disjoint_regions = section_geometry.check_geometry_disjoint(polygons)
	if disjoint_regions:
	warnings.warn(
	"\nThe section geometry contains disjoint regions which is invalid for warping analysis.\n"
	"Please revise your geometry to ensure there is connectivity between all regions.\n"
	"Please see https://sectionproperties.readthedocs.io/en/latest/rst/analysis.html#warping-analysis for more "
	"information."
	)

section-properties/sectionproperties/pre/geometry.py

Lines 2025 to 2038 in ca86307

	def check_geometry_disjoint(lop: List[Polygon]) -> bool:
	"""
	Returns True if any polygons in 'lop' are disjoint. Returns
	False, otherwise.
	"""
	# Build polygon connectivity network
	network = {}
	for idx_i, poly1 in enumerate(lop):
	for idx_j, poly2 in enumerate(lop):
	if idx_i != idx_j:
	connectivity = network.get(idx_i, set())
	if poly1.intersection(poly2):
	connectivity.add(idx_j)
	network[idx_i] = connectivity

Essentially it checks the intersection of several Shapely polygons. Generally, creating a polygon with zero area (and dealing with zero on computers in general) is complex, and can have unexpected results. Rather than setting the diameter to zero, I would try not creating those polygons at all.

[EngDM]

Oh that makes sense. Changing my code to check whether a diameter was specified as non-zero and only adding the reinforcing accordingly made the error resolve.

[EngDM]

When running the composite analysis am I correct in that to find the warping and torsion constants you would simply divide the EJ and EIw values by whatever your effective E is? In the documentation it 100% clear whether the E term is the effective E based on the symbol, although it does say weighted. I just want to make sure I'm not making a mistake.

[robbievanleeuwen]

For a composite analysis sectionproperties reports the torsional rigidity E.J. If you want to normalise this to a single material you can divide by that material's elastic modulus to obtain an effective J_eff.

Have a look at the concrete example above your screenshot.

[EngDM]

Yea that's what I assumed it would be. This whole joist reinforcing is a giant rabbit hole I swear. My governing code says to use section modulus to calculate moment, but for two materials with the same E but different Fy's it isn't very clear on how to calculate this. The only reliable way I've found is to just take the lower of the two Fy's and treat the composite shape as such. And then you have the whole issue of preloaded joists, and it's recommended to shore joists to return the stress in the existing member to near zero before installing reinforcing; but no contractor is going to do this for a small RTU job.

I can't believe what started out as just needing warping properties for a composite shape turned into such a mess..