Skip to content

Instantly share code, notes, and snippets.

@stephensmitchell
Last active December 28, 2025 15:55
Show Gist options
  • Select an option

  • Save stephensmitchell/bcdb0f8528b6ea0c61cfb134b8b7dc14 to your computer and use it in GitHub Desktop.

Select an option

Save stephensmitchell/bcdb0f8528b6ea0c61cfb134b8b7dc14 to your computer and use it in GitHub Desktop.
# created with Alibre Script Genie by Stephen S. Mitchell, https://github.com/stephensmitchell
from __future__ import division
import sys
import math
from AlibreScript import *
ScriptName = "Batch Area Moments"
ScriptVersion = "1.0"
DIALOG_WIDTH = 400
CURVE_SEGMENTS = 144
VERTEX_TOLERANCE = 1e-6
CIRCUMFERENCE_TOLERANCE = 0.02
ARC_LENGTH_TOLERANCE = 0.05
CURVATURE_THRESHOLD = 1.05
def vec_subtract(a, b):
return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]
def vec_add(a, b):
return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]
def vec_scale(v, s):
return [v[0]*s, v[1]*s, v[2]*s]
def vec_cross(a, b):
return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]
def vec_dot(a, b):
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
def vec_length(v):
return math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
def vec_normalize(v):
L = vec_length(v)
if L < 1e-12:
return [0.0, 0.0, 0.0]
return [v[0]/L, v[1]/L, v[2]/L]
def vec_distance(a, b):
return math.sqrt((b[0]-a[0])**2 + (b[1]-a[1])**2 + (b[2]-a[2])**2)
def vec_midpoint(a, b):
return [(a[0]+b[0])/2.0, (a[1]+b[1])/2.0, (a[2]+b[2])/2.0]
def vec_lerp(a, b, t):
return [a[0] + t*(b[0]-a[0]), a[1] + t*(b[1]-a[1]), a[2] + t*(b[2]-a[2])]
def safe_sqrt(val):
if val < 0:
return 0.0
return math.sqrt(val)
def safe_div(num, denom, default=float('inf')):
if abs(denom) < 1e-12:
return default
return num / denom
class CircleProperties:
def __init__(self, radius, center_x=0.0, center_y=0.0):
self.shape_type = "Circle"
self.radius = abs(radius)
self.n = 0
r = self.radius
A = math.pi * r * r
I = math.pi * r**4 / 4.0
self.area = A
self.centroid_x = center_x
self.centroid_y = center_y
self.Ix_centroid = I
self.Iy_centroid = I
self.Ixy_centroid = 0.0
self.J_centroid = 2.0 * I
self.Ix_origin = I + A * center_y**2
self.Iy_origin = I + A * center_x**2
self.Ixy_origin = A * center_x * center_y
self.I_max = I
self.I_min = I
self.theta_principal = 0.0
self.theta_principal_deg = 0.0
self.rx = r / 2.0
self.ry = r / 2.0
self.rp = r / math.sqrt(2.0)
self.c_top = r
self.c_bottom = r
self.c_right = r
self.c_left = r
S = safe_div(I, r, 0.0)
self.Sx_top = S
self.Sx_bottom = S
self.Sy_right = S
self.Sy_left = S
self.Sx_min = S
self.Sy_min = S
class AnnulusProperties:
def __init__(self, outer_radius, inner_radius, center_x=0.0, center_y=0.0):
self.shape_type = "Annulus"
R = abs(outer_radius)
r = abs(inner_radius)
if R < r:
R, r = r, R
self.outer_radius = R
self.inner_radius = r
self.n = 0
A = math.pi * (R**2 - r**2)
I = math.pi * (R**4 - r**4) / 4.0
self.area = A
self.centroid_x = center_x
self.centroid_y = center_y
self.Ix_centroid = I
self.Iy_centroid = I
self.Ixy_centroid = 0.0
self.J_centroid = 2.0 * I
self.Ix_origin = I + A * center_y**2
self.Iy_origin = I + A * center_x**2
self.Ixy_origin = A * center_x * center_y
self.I_max = I
self.I_min = I
self.theta_principal = 0.0
self.theta_principal_deg = 0.0
self.rx = safe_sqrt(safe_div(I, A, 0.0))
self.ry = safe_sqrt(safe_div(I, A, 0.0))
self.rp = safe_sqrt(safe_div(2.0 * I, A, 0.0))
self.c_top = R
self.c_bottom = R
self.c_right = R
self.c_left = R
S = safe_div(I, R, 0.0)
self.Sx_top = S
self.Sx_bottom = S
self.Sy_right = S
self.Sy_left = S
self.Sx_min = S
self.Sy_min = S
class RectangleProperties:
def __init__(self, width, height, center_x=0.0, center_y=0.0, rotation=0.0):
self.shape_type = "Rectangle"
self.width = abs(width)
self.height = abs(height)
self.rotation = rotation
self.n = 4
b = self.width
h = self.height
A = b * h
Ix_local = b * h**3 / 12.0
Iy_local = b**3 * h / 12.0
self.area = A
self.centroid_x = center_x
self.centroid_y = center_y
cos2 = math.cos(2.0 * rotation)
sin2 = math.sin(2.0 * rotation)
I_avg = (Ix_local + Iy_local) / 2.0
I_diff = (Ix_local - Iy_local) / 2.0
self.Ix_centroid = I_avg + I_diff * cos2
self.Iy_centroid = I_avg - I_diff * cos2
self.Ixy_centroid = -I_diff * sin2
self.J_centroid = Ix_local + Iy_local
self.Ix_origin = self.Ix_centroid + A * center_y**2
self.Iy_origin = self.Iy_centroid + A * center_x**2
self.Ixy_origin = self.Ixy_centroid + A * center_x * center_y
self.I_max = max(Ix_local, Iy_local)
self.I_min = min(Ix_local, Iy_local)
self.theta_principal = rotation if Ix_local >= Iy_local else rotation + math.pi/2.0
self.theta_principal_deg = math.degrees(self.theta_principal)
self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
cos_r = abs(math.cos(rotation))
sin_r = abs(math.sin(rotation))
self.c_top = (h * cos_r + b * sin_r) / 2.0
self.c_bottom = self.c_top
self.c_right = (b * cos_r + h * sin_r) / 2.0
self.c_left = self.c_right
self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
self.Sx_bottom = self.Sx_top
self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
self.Sy_left = self.Sy_right
self.Sx_min = self.Sx_top
self.Sy_min = self.Sy_right
class SemicircleProperties:
def __init__(self, radius, center_x=0.0, center_y=0.0, rotation=0.0):
self.shape_type = "Semicircle"
self.radius = abs(radius)
self.rotation = rotation
self.n = 0
r = self.radius
A = math.pi * r**2 / 2.0
y_bar = 4.0 * r / (3.0 * math.pi)
Ix_local = (math.pi/8.0 - 8.0/(9.0*math.pi)) * r**4
Iy_local = math.pi * r**4 / 8.0
self.area = A
cos_r = math.cos(rotation)
sin_r = math.sin(rotation)
self.centroid_x = center_x + y_bar * sin_r
self.centroid_y = center_y + y_bar * cos_r
cos2 = math.cos(2.0 * rotation)
sin2 = math.sin(2.0 * rotation)
I_avg = (Ix_local + Iy_local) / 2.0
I_diff = (Ix_local - Iy_local) / 2.0
self.Ix_centroid = I_avg + I_diff * cos2
self.Iy_centroid = I_avg - I_diff * cos2
self.Ixy_centroid = -I_diff * sin2
self.J_centroid = Ix_local + Iy_local
self.Ix_origin = self.Ix_centroid + A * self.centroid_y**2
self.Iy_origin = self.Iy_centroid + A * self.centroid_x**2
self.Ixy_origin = self.Ixy_centroid + A * self.centroid_x * self.centroid_y
self.I_max = max(Ix_local, Iy_local)
self.I_min = min(Ix_local, Iy_local)
self.theta_principal = rotation
self.theta_principal_deg = math.degrees(rotation)
self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
self.c_top = r - y_bar
self.c_bottom = y_bar
self.c_right = r
self.c_left = r
self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
self.Sx_min = min(self.Sx_top, self.Sx_bottom)
self.Sy_min = min(self.Sy_right, self.Sy_left)
class PolygonProperties:
def __init__(self, vertices_2d, shape_name="Polygon"):
self.shape_type = shape_name
self.vertices = vertices_2d
self.n = len(vertices_2d)
if self.n < 3:
raise ValueError("Polygon must have at least 3 vertices")
self._compute_all()
def _compute_all(self):
verts = self.vertices
n = self.n
signed_area = 0.0
for i in range(n):
j = (i + 1) % n
signed_area += verts[i][0] * verts[j][1]
signed_area -= verts[j][0] * verts[i][1]
signed_area *= 0.5
self.area = abs(signed_area)
if self.area < 1e-12:
raise ValueError("Polygon has zero or negligible area")
cx, cy = 0.0, 0.0
for i in range(n):
j = (i + 1) % n
cross = verts[i][0] * verts[j][1] - verts[j][0] * verts[i][1]
cx += (verts[i][0] + verts[j][0]) * cross
cy += (verts[i][1] + verts[j][1]) * cross
factor = 1.0 / (6.0 * signed_area) if abs(signed_area) > 1e-12 else 0.0
self.centroid_x = cx * factor
self.centroid_y = cy * factor
Ix, Iy, Ixy = 0.0, 0.0, 0.0
for i in range(n):
j = (i + 1) % n
x0, y0 = verts[i]
x1, y1 = verts[j]
cross = x0 * y1 - x1 * y0
Ix += (y0*y0 + y0*y1 + y1*y1) * cross
Iy += (x0*x0 + x0*x1 + x1*x1) * cross
Ixy += (x0*y1 + 2.0*x0*y0 + 2.0*x1*y1 + x1*y0) * cross
self.Ix_origin = abs(Ix / 12.0)
self.Iy_origin = abs(Iy / 12.0)
self.Ixy_origin = Ixy / 24.0 if signed_area > 0 else -Ixy / 24.0
A = self.area
self.Ix_centroid = abs(self.Ix_origin - A * self.centroid_y**2)
self.Iy_centroid = abs(self.Iy_origin - A * self.centroid_x**2)
self.Ixy_centroid = self.Ixy_origin - A * self.centroid_x * self.centroid_y
self.J_centroid = self.Ix_centroid + self.Iy_centroid
Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
I_avg = (Ix_c + Iy_c) / 2.0
I_diff = (Ix_c - Iy_c) / 2.0
R = safe_sqrt(I_diff**2 + Ixy_c**2)
self.I_max = I_avg + R
self.I_min = max(0.0, I_avg - R)
if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
self.theta_principal = 0.0
else:
self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
self.theta_principal_deg = math.degrees(self.theta_principal)
self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
x_rel = [v[0] - self.centroid_x for v in verts]
y_rel = [v[1] - self.centroid_y for v in verts]
self.c_top = max(y_rel) if y_rel else 0.0
self.c_bottom = abs(min(y_rel)) if y_rel else 0.0
self.c_right = max(x_rel) if x_rel else 0.0
self.c_left = abs(min(x_rel)) if x_rel else 0.0
self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
self.Sx_min = min(self.Sx_top, self.Sx_bottom)
self.Sy_min = min(self.Sy_right, self.Sy_left)
class PolygonWithCircularHoleProperties:
def __init__(self, outer_vertices_2d, hole_radius, hole_center_x=None, hole_center_y=None, shape_name="Polygon with Hole"):
self.shape_type = shape_name
outer = PolygonProperties(outer_vertices_2d, "outer")
if hole_center_x is None:
hole_center_x = outer.centroid_x
if hole_center_y is None:
hole_center_y = outer.centroid_y
r = abs(hole_radius)
hole_area = math.pi * r * r
hole_Ic = math.pi * r**4 / 4.0
dx = hole_center_x - outer.centroid_x
dy = hole_center_y - outer.centroid_y
hole_Ix_at_centroid = hole_Ic + hole_area * dy**2
hole_Iy_at_centroid = hole_Ic + hole_area * dx**2
hole_Ixy_at_centroid = hole_area * dx * dy
self.n = outer.n
self.area = outer.area - hole_area
if self.area < 1e-12:
raise ValueError("Hole area exceeds outer polygon area")
self.centroid_x = outer.centroid_x
self.centroid_y = outer.centroid_y
self.Ix_centroid = outer.Ix_centroid - hole_Ix_at_centroid
self.Iy_centroid = outer.Iy_centroid - hole_Iy_at_centroid
self.Ixy_centroid = outer.Ixy_centroid - hole_Ixy_at_centroid
self.J_centroid = self.Ix_centroid + self.Iy_centroid
self.Ix_origin = outer.Ix_origin - (hole_Ic + hole_area * hole_center_y**2)
self.Iy_origin = outer.Iy_origin - (hole_Ic + hole_area * hole_center_x**2)
self.Ixy_origin = outer.Ixy_origin - hole_area * hole_center_x * hole_center_y
Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
I_avg = (Ix_c + Iy_c) / 2.0
I_diff = (Ix_c - Iy_c) / 2.0
R = safe_sqrt(I_diff**2 + Ixy_c**2)
self.I_max = I_avg + R
self.I_min = max(0.0, I_avg - R)
if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
self.theta_principal = 0.0
else:
self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
self.theta_principal_deg = math.degrees(self.theta_principal)
self.rx = safe_sqrt(safe_div(self.Ix_centroid, self.area, 0.0))
self.ry = safe_sqrt(safe_div(self.Iy_centroid, self.area, 0.0))
self.rp = safe_sqrt(safe_div(self.J_centroid, self.area, 0.0))
self.c_top = outer.c_top
self.c_bottom = outer.c_bottom
self.c_right = outer.c_right
self.c_left = outer.c_left
self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
self.Sx_min = min(self.Sx_top, self.Sx_bottom)
self.Sy_min = min(self.Sy_right, self.Sy_left)
self.outer_area = outer.area
self.hole_radius = r
self.hole_area = hole_area
x_coords = [v[0] for v in outer_vertices_2d]
y_coords = [v[1] for v in outer_vertices_2d]
self.outer_width = max(x_coords) - min(x_coords)
self.outer_height = max(y_coords) - min(y_coords)
class EdgeInfo:
def __init__(self, edge, index):
self.edge = edge
self.index = index
self.diameter = None
self.length = None
self.vertices = []
self.vertex_count = 0
self.is_circular = False
self.is_closed = False
self.is_full_circle = False
self.is_curved = False
self.chord_length = 0.0
self._analyze()
def _analyze(self):
try:
d = self.edge.Diameter
if d is not None and d > 0:
self.diameter = d
self.is_circular = True
except Exception:
pass
try:
self.length = self.edge.Length
except Exception:
pass
try:
verts = self.edge.GetVertices()
if verts:
self.vertex_count = len(verts)
for v in verts:
self.vertices.append([v.X, v.Y, v.Z])
except Exception:
pass
self.is_closed = (self.vertex_count == 0)
if self.is_circular and self.diameter:
if self.is_closed:
self.is_full_circle = True
elif not self.length or self.length < 1e-9:
self.is_full_circle = True
elif self.length > 0:
expected_circumference = math.pi * self.diameter
if abs(self.length - expected_circumference) < expected_circumference * CIRCUMFERENCE_TOLERANCE:
self.is_full_circle = True
if self.vertex_count >= 2:
self.chord_length = vec_distance(self.vertices[0], self.vertices[1])
if self.length and self.chord_length > 1e-12:
ratio = self.length / self.chord_length
if ratio > CURVATURE_THRESHOLD:
self.is_curved = True
def is_complete_circle(edge_info):
return edge_info.is_closed or edge_info.is_full_circle
def collect_unique_vertices(edges, tolerance=VERTEX_TOLERANCE):
all_verts = []
for e in edges:
for v in e.vertices:
is_dup = False
for existing in all_verts:
if vec_distance(v, existing) < tolerance:
is_dup = True
break
if not is_dup:
all_verts.append(v)
return all_verts
def analyze_face_geometry(face, num_segments=CURVE_SEGMENTS):
alibre_area = None
try:
alibre_area = face.GetArea()
if alibre_area == 0:
alibre_area = None
except Exception:
pass
edges = []
try:
edges = face.GetEdges() or []
except Exception:
pass
face_vertices = []
try:
face_vertices = face.GetVertices() or []
except Exception:
pass
is_rect = False
try:
is_rect = face.IsRectangle()
except Exception:
pass
edge_infos = [EdgeInfo(e, i) for i, e in enumerate(edges)]
circular_edges = [e for e in edge_infos if e.is_circular]
curved_edges = [e for e in edge_infos if e.is_curved and not e.is_circular]
linear_edges = [e for e in edge_infos if not e.is_curved and not e.is_circular and e.vertex_count >= 2]
diameters = [e.diameter for e in circular_edges if e.diameter]
unique_diameters = list(set([round(d, 6) for d in diameters]))
if is_rect and len(linear_edges) == 4:
lengths = sorted([e.length for e in linear_edges if e.length])
if len(lengths) == 4:
width = (lengths[0] + lengths[1]) / 2.0
height = (lengths[2] + lengths[3]) / 2.0
return ("rectangle", RectangleProperties(width, height), alibre_area)
if len(circular_edges) == 1 and is_complete_circle(circular_edges[0]) and len(linear_edges) == 0:
radius = circular_edges[0].diameter / 2.0
return ("circle", CircleProperties(radius), alibre_area)
if len(circular_edges) == 2 and all(is_complete_circle(e) for e in circular_edges):
r1 = circular_edges[0].diameter / 2.0
r2 = circular_edges[1].diameter / 2.0
outer_r, inner_r = max(r1, r2), min(r1, r2)
if outer_r - inner_r < 0.001:
return ("circle", CircleProperties(outer_r), alibre_area)
return ("annulus", AnnulusProperties(outer_r, inner_r), alibre_area)
if len(unique_diameters) == 1 and len(circular_edges) >= 1 and len(linear_edges) == 0:
radius = unique_diameters[0] / 2.0
total_arc = sum(e.length for e in circular_edges if e.length)
expected_circumference = math.pi * unique_diameters[0]
if abs(total_arc - expected_circumference) < expected_circumference * 0.02:
return ("circle", CircleProperties(radius), alibre_area)
if len(unique_diameters) == 2:
r1, r2 = max(unique_diameters) / 2.0, min(unique_diameters) / 2.0
if r1 - r2 < 0.001:
return ("circle", CircleProperties(r1), alibre_area)
return ("annulus", AnnulusProperties(r1, r2), alibre_area)
if len(circular_edges) == 1 and len(linear_edges) == 1:
arc_edge = circular_edges[0]
if arc_edge.diameter and arc_edge.length:
expected_semicircle = math.pi * arc_edge.diameter / 2.0
if abs(arc_edge.length - expected_semicircle) < expected_semicircle * ARC_LENGTH_TOLERANCE:
radius = arc_edge.diameter / 2.0
return ("semicircle", SemicircleProperties(radius), alibre_area)
if len(linear_edges) >= 3 and len(circular_edges) == 1 and is_complete_circle(circular_edges[0]):
hole_edge = circular_edges[0]
if hole_edge.diameter:
hole_radius = hole_edge.diameter / 2.0
outer_verts = collect_unique_vertices(linear_edges)
if len(outer_verts) >= 3:
outer_2d = project_to_2d(outer_verts)
shape_name = "Polygon (%d sides) with Hole" % len(linear_edges)
return ("polygon_with_hole", PolygonWithCircularHoleProperties(outer_2d, hole_radius, shape_name=shape_name), alibre_area)
if len(linear_edges) == 3 and len(circular_edges) == 0:
all_verts = collect_unique_vertices(linear_edges)
if len(all_verts) == 3:
points_2d = project_to_2d(all_verts)
return ("triangle", PolygonProperties(points_2d, "Triangle"), alibre_area)
if len(linear_edges) >= 3 and len(circular_edges) == 0 and len(curved_edges) == 0:
all_verts = collect_unique_vertices(linear_edges)
if len(all_verts) >= 3:
points_2d = project_to_2d(all_verts)
return ("polygon", PolygonProperties(points_2d, "Polygon (%d sides)" % len(linear_edges)), alibre_area)
if len(curved_edges) > 0 or (len(circular_edges) > 0 and not all(is_complete_circle(e) for e in circular_edges)):
boundary_points = discretize_face_boundary(edge_infos, num_segments)
if len(boundary_points) >= 3:
points_2d = project_to_2d(boundary_points)
return ("freeform", PolygonProperties(points_2d, "Freeform (%d pts)" % len(points_2d)), alibre_area)
boundary_points = []
for v in face_vertices:
try:
pt = [v.X, v.Y, v.Z]
is_dup = any(vec_distance(pt, ex) < VERTEX_TOLERANCE for ex in boundary_points)
if not is_dup:
boundary_points.append(pt)
except Exception:
pass
if len(boundary_points) >= 3:
points_2d = project_to_2d(boundary_points)
return ("polygon", PolygonProperties(points_2d, "Polygon (from vertices)"), alibre_area)
if alibre_area and alibre_area > 0:
radius = math.sqrt(alibre_area / math.pi)
props = CircleProperties(radius)
props.shape_type = "Equivalent Circle (from area)"
return ("equivalent_circle", props, alibre_area)
return ("unknown", None, alibre_area)
def discretize_face_boundary(edge_infos, num_segments):
boundary_points = []
ordered_edges = order_edges_by_connectivity(edge_infos)
for edge_info in ordered_edges:
if is_complete_circle(edge_info) and edge_info.diameter:
radius = edge_info.diameter / 2.0
center = estimate_circle_center(edge_info)
for i in range(num_segments):
angle = 2.0 * math.pi * i / num_segments
pt = [center[0] + radius * math.cos(angle), center[1] + radius * math.sin(angle), center[2]]
boundary_points.append(pt)
elif edge_info.is_curved or edge_info.is_circular:
pts = discretize_curved_edge(edge_info, max(8, num_segments // 4))
for pt in pts:
is_dup = any(vec_distance(pt, ex) < VERTEX_TOLERANCE for ex in boundary_points)
if not is_dup:
boundary_points.append(pt)
else:
for v in edge_info.vertices:
is_dup = any(vec_distance(v, ex) < VERTEX_TOLERANCE for ex in boundary_points)
if not is_dup:
boundary_points.append(v)
return boundary_points
def order_edges_by_connectivity(edge_infos):
if not edge_infos:
return []
ordered = []
remaining = list(edge_infos)
ordered.append(remaining.pop(0))
while remaining:
last_edge = ordered[-1]
found_idx = -1
for i, edge in enumerate(remaining):
if edges_share_vertex(last_edge, edge):
found_idx = i
break
if found_idx >= 0:
ordered.append(remaining.pop(found_idx))
else:
ordered.append(remaining.pop(0))
return ordered
def edges_share_vertex(edge1, edge2):
for v1 in edge1.vertices:
for v2 in edge2.vertices:
if vec_distance(v1, v2) < VERTEX_TOLERANCE:
return True
return False
def estimate_circle_center(edge_info):
if edge_info.vertices:
n = len(edge_info.vertices)
cx = sum(v[0] for v in edge_info.vertices) / n
cy = sum(v[1] for v in edge_info.vertices) / n
cz = sum(v[2] for v in edge_info.vertices) / n
return [cx, cy, cz]
return [0.0, 0.0, 0.0]
def discretize_curved_edge(edge_info, num_points):
points = []
num_points = max(4, num_points)
if len(edge_info.vertices) >= 2:
v1, v2 = edge_info.vertices[0], edge_info.vertices[1]
for i in range(num_points + 1):
t = i / float(num_points)
points.append(vec_lerp(v1, v2, t))
return points
def project_to_2d(points_3d):
n = len(points_3d)
if n < 3:
raise ValueError("Need at least 3 points for 2D projection")
cx = sum(p[0] for p in points_3d) / n
cy = sum(p[1] for p in points_3d) / n
cz = sum(p[2] for p in points_3d) / n
origin = [cx, cy, cz]
p0 = points_3d[0]
v1 = None
for i in range(1, n):
vec = vec_subtract(points_3d[i], p0)
if vec_length(vec) > VERTEX_TOLERANCE:
v1 = vec
break
if v1 is None:
raise ValueError("All points are coincident")
v2 = None
for i in range(2, n):
vec = vec_subtract(points_3d[i], p0)
cross = vec_cross(v1, vec)
if vec_length(cross) > VERTEX_TOLERANCE:
v2 = vec
break
if v2 is None:
if abs(v1[2]) < 0.9:
v2 = vec_cross(v1, [0.0, 0.0, 1.0])
else:
v2 = vec_cross(v1, [1.0, 0.0, 0.0])
normal = vec_normalize(vec_cross(v1, v2))
u_axis = vec_normalize(v1)
v_axis = vec_cross(normal, u_axis)
points_2d = []
for p3d in points_3d:
rel = vec_subtract(p3d, origin)
u = vec_dot(rel, u_axis)
v = vec_dot(rel, v_axis)
points_2d.append([u, v])
cx2 = sum(p[0] for p in points_2d) / n
cy2 = sum(p[1] for p in points_2d) / n
def angle_key(p):
return math.atan2(p[1] - cy2, p[0] - cx2)
points_2d.sort(key=angle_key)
return points_2d
def fmt(value, decimals=6):
if value is None:
return "N/A"
if isinstance(value, float) and (math.isinf(value) or math.isnan(value)):
return "N/A"
if abs(value) < 1e-10:
return "0.0"
elif abs(value) >= 1e6 or (abs(value) < 0.001 and abs(value) > 1e-10):
return "%.4e" % value
else:
return ("%%.%df" % decimals) % value
UNIT_CONVERSIONS = {
'mm': 1.0,
'cm': 0.1,
'm': 0.001,
'in': 1.0 / 25.4,
'ft': 1.0 / 304.8
}
def convert_value(value, power, unit):
if value is None or (isinstance(value, float) and (math.isinf(value) or math.isnan(value))):
return value
factor = UNIT_CONVERSIONS[unit] ** power
return value * factor
def generate_face_section(face_name, props, alibre_area, unit):
u1 = unit
u2 = unit + "^2"
u3 = unit + "^3"
u4 = unit + "^4"
p = props
lines = []
lines.append(" %s (%s)" % (face_name, props.shape_type))
lines.append(" " + "-" * 40)
if hasattr(props, 'radius') and not hasattr(props, 'outer_radius'):
lines.append(" Radius: %s %s" % (fmt(convert_value(props.radius, 1, unit)), u1))
if hasattr(props, 'outer_radius'):
lines.append(" Outer R: %s %s" % (fmt(convert_value(props.outer_radius, 1, unit)), u1))
lines.append(" Inner R: %s %s" % (fmt(convert_value(props.inner_radius, 1, unit)), u1))
if hasattr(props, 'width') and hasattr(props, 'height'):
lines.append(" Width: %s %s" % (fmt(convert_value(props.width, 1, unit)), u1))
lines.append(" Height: %s %s" % (fmt(convert_value(props.height, 1, unit)), u1))
if hasattr(props, 'outer_width') and hasattr(props, 'outer_height'):
lines.append(" Outer W: %s %s" % (fmt(convert_value(props.outer_width, 1, unit)), u1))
lines.append(" Outer H: %s %s" % (fmt(convert_value(props.outer_height, 1, unit)), u1))
if hasattr(props, 'hole_radius'):
lines.append(" Hole R: %s %s" % (fmt(convert_value(props.hole_radius, 1, unit)), u1))
lines.append(" Area: %s %s" % (fmt(convert_value(p.area, 2, unit)), u2))
lines.append(" Ix-x: %s %s" % (fmt(convert_value(p.Ix_centroid, 4, unit)), u4))
lines.append(" Iy-y: %s %s" % (fmt(convert_value(p.Iy_centroid, 4, unit)), u4))
lines.append(" J: %s %s" % (fmt(convert_value(p.J_centroid, 4, unit)), u4))
lines.append(" rx-x: %s %s" % (fmt(convert_value(p.rx, 1, unit)), u1))
lines.append(" ry-y: %s %s" % (fmt(convert_value(p.ry, 1, unit)), u1))
lines.append(" Sx-x (min): %s %s" % (fmt(convert_value(p.Sx_min, 3, unit)), u3))
lines.append(" Sy-y (min): %s %s" % (fmt(convert_value(p.Sy_min, 3, unit)), u3))
lines.append("")
return "\n".join(lines)
def generate_batch_report(face_results):
lines = []
lines.append("")
lines.append("=" * 70)
lines.append(" BATCH AREA MOMENTS REPORT")
lines.append(" %s v%s" % (ScriptName, ScriptVersion))
lines.append("=" * 70)
lines.append("")
lines.append("Faces Analyzed: %d" % len(face_results))
lines.append("")
for unit in ['mm', 'cm', 'm', 'in', 'ft']:
lines.append("=" * 70)
lines.append("RESULTS IN %s" % unit.upper())
lines.append("=" * 70)
lines.append("")
for face_name, props, alibre_area in face_results:
if props is not None:
lines.append(generate_face_section(face_name, props, alibre_area, unit))
lines.append("=" * 70)
lines.append("END OF REPORT")
lines.append("=" * 70)
return "\n".join(lines)
def get_all_faces(part):
faces = []
try:
all_faces = part.GetFaces()
if all_faces:
for face in all_faces:
try:
name = face.Name if hasattr(face, 'Name') and face.Name else None
if name:
faces.append((name, face))
except Exception:
pass
except Exception:
pass
faces.sort(key=lambda x: x[0])
return faces
def run():
Win = Windows()
try:
part = CurrentPart()
if part is None:
Win.ErrorDialog("Please open a part before running this script.", ScriptName)
return
except Exception:
Win.ErrorDialog("Please open a part before running this script.", ScriptName)
return
all_faces = get_all_faces(part)
if not all_faces:
Win.ErrorDialog("No named faces found in the part.\n\nFaces should be named like Face<1>, Face<2>, etc.", ScriptName)
return
print "=" * 70
print "BATCH AREA MOMENTS - Scanning %d faces..." % len(all_faces)
print "=" * 70
face_results = []
for face_name, face in all_faces:
print "Analyzing: %s" % face_name
try:
shape_type, props, alibre_area = analyze_face_geometry(face, CURVE_SEGMENTS)
if props is not None:
if shape_type not in ("polygon_with_hole", "circle", "annulus", "rectangle"):
if alibre_area and alibre_area > 0 and props.area > 0:
area_ratio = alibre_area / props.area
if area_ratio > 1.10 or area_ratio < 0.90:
print " -> Skipping (non-planar, area mismatch: %.1f%%)" % ((area_ratio - 1.0) * 100)
continue
face_results.append((face_name, props, alibre_area))
print " -> %s (Area: %.2f mm^2)" % (props.shape_type, props.area)
else:
print " -> Could not determine geometry"
except Exception as ex:
print " -> Error: %s" % str(ex)
print ""
if not face_results:
Win.ErrorDialog("No faces could be analyzed.\n\nMake sure faces are planar.", ScriptName)
return
report = generate_batch_report(face_results)
print report
sys.stdout.flush()
Win.InfoDialog("Batch analysis complete!\n\n%d faces analyzed.\nSee console for full report." % len(face_results), ScriptName)
run()
Topology Summary:
Type Total Failed
------------ ----- ------
Bodies 1 0
Faces 10 0
Edges 16 0
Vertices 12 0
Lumps 1 -
Shells 1 -
Loops 14 -
Coedges 32 -
TEdges 0 -
TVertices 0 -
TCoedges 0 -
======================================================================
BATCH AREA MOMENTS - Scanning 10 faces...
======================================================================
Analyzing: Face<10>
-> Polygon (4 sides) (Area: 13664.78 mm^2)
Analyzing: Face<1>
-> Annulus (Area: 335.64 mm^2)
Analyzing: Face<2>
-> Circle (Area: 1140.09 mm^2)
Analyzing: Face<3>
-> Circle (Area: 804.45 mm^2)
Analyzing: Face<4>
-> Annulus (Area: 2026.83 mm^2)
Analyzing: Face<5>
-> Polygon (4 sides) with Hole (Area: 2639.52 mm^2)
Analyzing: Face<6>
-> Rectangle (Area: 47096.73 mm^2)
Analyzing: Face<7>
-> Polygon (4 sides) (Area: 13643.50 mm^2)
Analyzing: Face<8>
-> Polygon (4 sides) (Area: 13664.78 mm^2)
Analyzing: Face<9>
-> Polygon (4 sides) (Area: 13643.50 mm^2)
======================================================================
BATCH AREA MOMENTS REPORT
Batch Area Moments v1.0
======================================================================
Faces Analyzed: 10
======================================================================
RESULTS IN MM
======================================================================
Face<10> (Polygon (4 sides))
----------------------------------------
Area: 13664.781537 mm^2
Ix-x: 9.2139e+06 mm^4
Iy-y: 2.9826e+07 mm^4
J: 3.9040e+07 mm^4
rx-x: 25.966873 mm
ry-y: 46.719223 mm
Sx-x (min): 169835.357191 mm^3
Sy-y (min): 276397.051143 mm^3
Face<1> (Annulus)
----------------------------------------
Outer R: 19.050000 mm
Inner R: 16.002000 mm
Area: 335.643034 mm^2
Ix-x: 51937.948861 mm^4
Iy-y: 51937.948861 mm^4
J: 103875.897721 mm^4
rx-x: 12.439519 mm
ry-y: 12.439519 mm
Sx-x (min): 2726.401515 mm^3
Sy-y (min): 2726.401515 mm^3
Face<2> (Circle)
----------------------------------------
Radius: 19.050000 mm
Area: 1140.091828 mm^2
Ix-x: 103435.543650 mm^4
Iy-y: 103435.543650 mm^4
J: 206871.087300 mm^4
rx-x: 9.525000 mm
ry-y: 9.525000 mm
Sx-x (min): 5429.687331 mm^3
Sy-y (min): 5429.687331 mm^3
Face<3> (Circle)
----------------------------------------
Radius: 16.002000 mm
Area: 804.448794 mm^2
Ix-x: 51497.594789 mm^4
Iy-y: 51497.594789 mm^4
J: 102995.189579 mm^4
rx-x: 8.001000 mm
ry-y: 8.001000 mm
Sx-x (min): 3218.197400 mm^3
Sy-y (min): 3218.197400 mm^3
Face<4> (Annulus)
----------------------------------------
Outer R: 31.750000 mm
Inner R: 19.050000 mm
Area: 2026.829916 mm^2
Ix-x: 694678.219081 mm^4
Iy-y: 694678.219081 mm^4
J: 1.3894e+06 mm^4
rx-x: 18.513272 mm
ry-y: 18.513272 mm
Sx-x (min): 21879.628947 mm^3
Sy-y (min): 21879.628947 mm^3
Face<5> (Polygon (4 sides) with Hole)
----------------------------------------
Outer W: 76.200000 mm
Outer H: 76.200000 mm
Hole R: 31.750000 mm
Area: 2639.518256 mm^2
Ix-x: 2.0114e+06 mm^4
Iy-y: 2.0114e+06 mm^4
J: 4.0229e+06 mm^4
rx-x: 27.605277 mm
ry-y: 27.605277 mm
Sx-x (min): 52793.920212 mm^3
Sy-y (min): 52793.920212 mm^3
Face<6> (Rectangle)
----------------------------------------
Width: 215.819414 mm
Height: 218.222875 mm
Area: 47096.733130 mm^2
Ix-x: 1.8690e+08 mm^4
Iy-y: 1.8281e+08 mm^4
J: 3.6971e+08 mm^4
rx-x: 62.995518 mm
ry-y: 62.301698 mm
Sx-x (min): 1.7129e+06 mm^3
Sy-y (min): 1.6941e+06 mm^3
Face<7> (Polygon (4 sides))
----------------------------------------
Area: 13643.503951 mm^2
Ix-x: 9.0085e+06 mm^4
Iy-y: 3.0373e+07 mm^4
J: 3.9381e+07 mm^4
rx-x: 25.695803 mm
ry-y: 47.182122 mm
Sx-x (min): 167471.687718 mm^3
Sy-y (min): 278362.405971 mm^3
Face<8> (Polygon (4 sides))
----------------------------------------
Area: 13664.781537 mm^2
Ix-x: 9.2139e+06 mm^4
Iy-y: 2.9826e+07 mm^4
J: 3.9040e+07 mm^4
rx-x: 25.966873 mm
ry-y: 46.719223 mm
Sx-x (min): 169835.357191 mm^3
Sy-y (min): 276397.051143 mm^3
Face<9> (Polygon (4 sides))
----------------------------------------
Area: 13643.503951 mm^2
Ix-x: 9.0085e+06 mm^4
Iy-y: 3.0373e+07 mm^4
J: 3.9381e+07 mm^4
rx-x: 25.695803 mm
ry-y: 47.182122 mm
Sx-x (min): 167471.687718 mm^3
Sy-y (min): 278362.405971 mm^3
======================================================================
RESULTS IN CM
======================================================================
Face<10> (Polygon (4 sides))
----------------------------------------
Area: 136.647815 cm^2
Ix-x: 921.386819 cm^4
Iy-y: 2982.592485 cm^4
J: 3903.979304 cm^4
rx-x: 2.596687 cm
ry-y: 4.671922 cm
Sx-x (min): 169.835357 cm^3
Sy-y (min): 276.397051 cm^3
Face<1> (Annulus)
----------------------------------------
Outer R: 1.905000 cm
Inner R: 1.600200 cm
Area: 3.356430 cm^2
Ix-x: 5.193795 cm^4
Iy-y: 5.193795 cm^4
J: 10.387590 cm^4
rx-x: 1.243952 cm
ry-y: 1.243952 cm
Sx-x (min): 2.726402 cm^3
Sy-y (min): 2.726402 cm^3
Face<2> (Circle)
----------------------------------------
Radius: 1.905000 cm
Area: 11.400918 cm^2
Ix-x: 10.343554 cm^4
Iy-y: 10.343554 cm^4
J: 20.687109 cm^4
rx-x: 0.952500 cm
ry-y: 0.952500 cm
Sx-x (min): 5.429687 cm^3
Sy-y (min): 5.429687 cm^3
Face<3> (Circle)
----------------------------------------
Radius: 1.600200 cm
Area: 8.044488 cm^2
Ix-x: 5.149759 cm^4
Iy-y: 5.149759 cm^4
J: 10.299519 cm^4
rx-x: 0.800100 cm
ry-y: 0.800100 cm
Sx-x (min): 3.218197 cm^3
Sy-y (min): 3.218197 cm^3
Face<4> (Annulus)
----------------------------------------
Outer R: 3.175000 cm
Inner R: 1.905000 cm
Area: 20.268299 cm^2
Ix-x: 69.467822 cm^4
Iy-y: 69.467822 cm^4
J: 138.935644 cm^4
rx-x: 1.851327 cm
ry-y: 1.851327 cm
Sx-x (min): 21.879629 cm^3
Sy-y (min): 21.879629 cm^3
Face<5> (Polygon (4 sides) with Hole)
----------------------------------------
Outer W: 7.620000 cm
Outer H: 7.620000 cm
Hole R: 3.175000 cm
Area: 26.395183 cm^2
Ix-x: 201.144836 cm^4
Iy-y: 201.144836 cm^4
J: 402.289672 cm^4
rx-x: 2.760528 cm
ry-y: 2.760528 cm
Sx-x (min): 52.793920 cm^3
Sy-y (min): 52.793920 cm^3
Face<6> (Rectangle)
----------------------------------------
Width: 21.581941 cm
Height: 21.822288 cm
Area: 470.967331 cm^2
Ix-x: 18690.033702 cm^4
Iy-y: 18280.604657 cm^4
J: 36970.638359 cm^4
rx-x: 6.299552 cm
ry-y: 6.230170 cm
Sx-x (min): 1712.930753 cm^3
Sy-y (min): 1694.064893 cm^3
Face<7> (Polygon (4 sides))
----------------------------------------
Area: 136.435040 cm^2
Ix-x: 900.845517 cm^4
Iy-y: 3037.252230 cm^4
J: 3938.097746 cm^4
rx-x: 2.569580 cm
ry-y: 4.718212 cm
Sx-x (min): 167.471688 cm^3
Sy-y (min): 278.362406 cm^3
Face<8> (Polygon (4 sides))
----------------------------------------
Area: 136.647815 cm^2
Ix-x: 921.386819 cm^4
Iy-y: 2982.592485 cm^4
J: 3903.979304 cm^4
rx-x: 2.596687 cm
ry-y: 4.671922 cm
Sx-x (min): 169.835357 cm^3
Sy-y (min): 276.397051 cm^3
Face<9> (Polygon (4 sides))
----------------------------------------
Area: 136.435040 cm^2
Ix-x: 900.845517 cm^4
Iy-y: 3037.252230 cm^4
J: 3938.097746 cm^4
rx-x: 2.569580 cm
ry-y: 4.718212 cm
Sx-x (min): 167.471688 cm^3
Sy-y (min): 278.362406 cm^3
======================================================================
RESULTS IN M
======================================================================
Face<10> (Polygon (4 sides))
----------------------------------------
Area: 0.013665 m^2
Ix-x: 9.2139e-06 m^4
Iy-y: 2.9826e-05 m^4
J: 3.9040e-05 m^4
rx-x: 0.025967 m
ry-y: 0.046719 m
Sx-x (min): 1.6984e-04 m^3
Sy-y (min): 2.7640e-04 m^3
Face<1> (Annulus)
----------------------------------------
Outer R: 0.019050 m
Inner R: 0.016002 m
Area: 3.3564e-04 m^2
Ix-x: 5.1938e-08 m^4
Iy-y: 5.1938e-08 m^4
J: 1.0388e-07 m^4
rx-x: 0.012440 m
ry-y: 0.012440 m
Sx-x (min): 2.7264e-06 m^3
Sy-y (min): 2.7264e-06 m^3
Face<2> (Circle)
----------------------------------------
Radius: 0.019050 m
Area: 0.001140 m^2
Ix-x: 1.0344e-07 m^4
Iy-y: 1.0344e-07 m^4
J: 2.0687e-07 m^4
rx-x: 0.009525 m
ry-y: 0.009525 m
Sx-x (min): 5.4297e-06 m^3
Sy-y (min): 5.4297e-06 m^3
Face<3> (Circle)
----------------------------------------
Radius: 0.016002 m
Area: 8.0445e-04 m^2
Ix-x: 5.1498e-08 m^4
Iy-y: 5.1498e-08 m^4
J: 1.0300e-07 m^4
rx-x: 0.008001 m
ry-y: 0.008001 m
Sx-x (min): 3.2182e-06 m^3
Sy-y (min): 3.2182e-06 m^3
Face<4> (Annulus)
----------------------------------------
Outer R: 0.031750 m
Inner R: 0.019050 m
Area: 0.002027 m^2
Ix-x: 6.9468e-07 m^4
Iy-y: 6.9468e-07 m^4
J: 1.3894e-06 m^4
rx-x: 0.018513 m
ry-y: 0.018513 m
Sx-x (min): 2.1880e-05 m^3
Sy-y (min): 2.1880e-05 m^3
Face<5> (Polygon (4 sides) with Hole)
----------------------------------------
Outer W: 0.076200 m
Outer H: 0.076200 m
Hole R: 0.031750 m
Area: 0.002640 m^2
Ix-x: 2.0114e-06 m^4
Iy-y: 2.0114e-06 m^4
J: 4.0229e-06 m^4
rx-x: 0.027605 m
ry-y: 0.027605 m
Sx-x (min): 5.2794e-05 m^3
Sy-y (min): 5.2794e-05 m^3
Face<6> (Rectangle)
----------------------------------------
Width: 0.215819 m
Height: 0.218223 m
Area: 0.047097 m^2
Ix-x: 1.8690e-04 m^4
Iy-y: 1.8281e-04 m^4
J: 3.6971e-04 m^4
rx-x: 0.062996 m
ry-y: 0.062302 m
Sx-x (min): 0.001713 m^3
Sy-y (min): 0.001694 m^3
Face<7> (Polygon (4 sides))
----------------------------------------
Area: 0.013644 m^2
Ix-x: 9.0085e-06 m^4
Iy-y: 3.0373e-05 m^4
J: 3.9381e-05 m^4
rx-x: 0.025696 m
ry-y: 0.047182 m
Sx-x (min): 1.6747e-04 m^3
Sy-y (min): 2.7836e-04 m^3
Face<8> (Polygon (4 sides))
----------------------------------------
Area: 0.013665 m^2
Ix-x: 9.2139e-06 m^4
Iy-y: 2.9826e-05 m^4
J: 3.9040e-05 m^4
rx-x: 0.025967 m
ry-y: 0.046719 m
Sx-x (min): 1.6984e-04 m^3
Sy-y (min): 2.7640e-04 m^3
Face<9> (Polygon (4 sides))
----------------------------------------
Area: 0.013644 m^2
Ix-x: 9.0085e-06 m^4
Iy-y: 3.0373e-05 m^4
J: 3.9381e-05 m^4
rx-x: 0.025696 m
ry-y: 0.047182 m
Sx-x (min): 1.6747e-04 m^3
Sy-y (min): 2.7836e-04 m^3
======================================================================
RESULTS IN IN
======================================================================
Face<10> (Polygon (4 sides))
----------------------------------------
Area: 21.180454 in^2
Ix-x: 22.136407 in^4
Iy-y: 71.657071 in^4
J: 93.793478 in^4
rx-x: 1.022318 in
ry-y: 1.839339 in
Sx-x (min): 10.363989 in^3
Sy-y (min): 16.866783 in^3
Face<1> (Annulus)
----------------------------------------
Outer R: 0.750000 in
Inner R: 0.630000 in
Area: 0.520248 in^2
Ix-x: 0.124781 in^4
Iy-y: 0.124781 in^4
J: 0.249563 in^4
rx-x: 0.489745 in
ry-y: 0.489745 in
Sx-x (min): 0.166375 in^3
Sy-y (min): 0.166375 in^3
Face<2> (Circle)
----------------------------------------
Radius: 0.750000 in
Area: 1.767146 in^2
Ix-x: 0.248505 in^4
Iy-y: 0.248505 in^4
J: 0.497010 in^4
rx-x: 0.375000 in
ry-y: 0.375000 in
Sx-x (min): 0.331340 in^3
Sy-y (min): 0.331340 in^3
Face<3> (Circle)
----------------------------------------
Radius: 0.630000 in
Area: 1.246898 in^2
Ix-x: 0.123723 in^4
Iy-y: 0.123723 in^4
J: 0.247447 in^4
rx-x: 0.315000 in
ry-y: 0.315000 in
Sx-x (min): 0.196386 in^3
Sy-y (min): 0.196386 in^3
Face<4> (Annulus)
----------------------------------------
Outer R: 1.250000 in
Inner R: 0.750000 in
Area: 3.141593 in^2
Ix-x: 1.668971 in^4
Iy-y: 1.668971 in^4
J: 3.337942 in^4
rx-x: 0.728869 in
ry-y: 0.728869 in
Sx-x (min): 1.335177 in^3
Sy-y (min): 1.335177 in^3
Face<5> (Polygon (4 sides) with Hole)
----------------------------------------
Outer W: 3.000000 in
Outer H: 3.000000 in
Hole R: 1.250000 in
Area: 4.091261 in^2
Ix-x: 4.832524 in^4
Iy-y: 4.832524 in^4
J: 9.665048 in^4
rx-x: 1.086822 in
ry-y: 1.086822 in
Sx-x (min): 3.221683 in^3
Sy-y (min): 3.221683 in^3
Face<6> (Rectangle)
----------------------------------------
Width: 8.496827 in
Height: 8.591452 in
Area: 73.000082 in^2
Ix-x: 449.029856 in^4
Iy-y: 439.193284 in^4
J: 888.223139 in^4
rx-x: 2.480138 in
ry-y: 2.452823 in
Sx-x (min): 104.529448 in^3
Sy-y (min): 103.378183 in^3
Face<7> (Polygon (4 sides))
----------------------------------------
Area: 21.147473 in^2
Ix-x: 21.642900 in^4
Iy-y: 72.970277 in^4
J: 94.613177 in^4
rx-x: 1.011646 in
ry-y: 1.857564 in
Sx-x (min): 10.219749 in^3
Sy-y (min): 16.986716 in^3
Face<8> (Polygon (4 sides))
----------------------------------------
Area: 21.180454 in^2
Ix-x: 22.136407 in^4
Iy-y: 71.657071 in^4
J: 93.793478 in^4
rx-x: 1.022318 in
ry-y: 1.839339 in
Sx-x (min): 10.363989 in^3
Sy-y (min): 16.866783 in^3
Face<9> (Polygon (4 sides))
----------------------------------------
Area: 21.147473 in^2
Ix-x: 21.642900 in^4
Iy-y: 72.970277 in^4
J: 94.613177 in^4
rx-x: 1.011646 in
ry-y: 1.857564 in
Sx-x (min): 10.219749 in^3
Sy-y (min): 16.986716 in^3
======================================================================
RESULTS IN FT
======================================================================
Face<10> (Polygon (4 sides))
----------------------------------------
Area: 0.147086 ft^2
Ix-x: 0.001068 ft^4
Iy-y: 0.003456 ft^4
J: 0.004523 ft^4
rx-x: 0.085193 ft
ry-y: 0.153278 ft
Sx-x (min): 0.005998 ft^3
Sy-y (min): 0.009761 ft^3
Face<1> (Annulus)
----------------------------------------
Outer R: 0.062500 ft
Inner R: 0.052500 ft
Area: 0.003613 ft^2
Ix-x: 6.0176e-06 ft^4
Iy-y: 6.0176e-06 ft^4
J: 1.2035e-05 ft^4
rx-x: 0.040812 ft
ry-y: 0.040812 ft
Sx-x (min): 9.6282e-05 ft^3
Sy-y (min): 9.6282e-05 ft^3
Face<2> (Circle)
----------------------------------------
Radius: 0.062500 ft
Area: 0.012272 ft^2
Ix-x: 1.1984e-05 ft^4
Iy-y: 1.1984e-05 ft^4
J: 2.3968e-05 ft^4
rx-x: 0.031250 ft
ry-y: 0.031250 ft
Sx-x (min): 1.9175e-04 ft^3
Sy-y (min): 1.9175e-04 ft^3
Face<3> (Circle)
----------------------------------------
Radius: 0.052500 ft
Area: 0.008659 ft^2
Ix-x: 5.9666e-06 ft^4
Iy-y: 5.9666e-06 ft^4
J: 1.1933e-05 ft^4
rx-x: 0.026250 ft
ry-y: 0.026250 ft
Sx-x (min): 1.1365e-04 ft^3
Sy-y (min): 1.1365e-04 ft^3
Face<4> (Annulus)
----------------------------------------
Outer R: 0.104167 ft
Inner R: 0.062500 ft
Area: 0.021817 ft^2
Ix-x: 8.0487e-05 ft^4
Iy-y: 8.0487e-05 ft^4
J: 1.6097e-04 ft^4
rx-x: 0.060739 ft
ry-y: 0.060739 ft
Sx-x (min): 7.7267e-04 ft^3
Sy-y (min): 7.7267e-04 ft^3
Face<5> (Polygon (4 sides) with Hole)
----------------------------------------
Outer W: 0.250000 ft
Outer H: 0.250000 ft
Hole R: 0.104167 ft
Area: 0.028412 ft^2
Ix-x: 2.3305e-04 ft^4
Iy-y: 2.3305e-04 ft^4
J: 4.6610e-04 ft^4
rx-x: 0.090568 ft
ry-y: 0.090568 ft
Sx-x (min): 0.001864 ft^3
Sy-y (min): 0.001864 ft^3
Face<6> (Rectangle)
----------------------------------------
Width: 0.708069 ft
Height: 0.715954 ft
Area: 0.506945 ft^2
Ix-x: 0.021655 ft^4
Iy-y: 0.021180 ft^4
J: 0.042835 ft^4
rx-x: 0.206678 ft
ry-y: 0.204402 ft
Sx-x (min): 0.060492 ft^3
Sy-y (min): 0.059825 ft^3
Face<7> (Polygon (4 sides))
----------------------------------------
Area: 0.146857 ft^2
Ix-x: 0.001044 ft^4
Iy-y: 0.003519 ft^4
J: 0.004563 ft^4
rx-x: 0.084304 ft
ry-y: 0.154797 ft
Sx-x (min): 0.005914 ft^3
Sy-y (min): 0.009830 ft^3
Face<8> (Polygon (4 sides))
----------------------------------------
Area: 0.147086 ft^2
Ix-x: 0.001068 ft^4
Iy-y: 0.003456 ft^4
J: 0.004523 ft^4
rx-x: 0.085193 ft
ry-y: 0.153278 ft
Sx-x (min): 0.005998 ft^3
Sy-y (min): 0.009761 ft^3
Face<9> (Polygon (4 sides))
----------------------------------------
Area: 0.146857 ft^2
Ix-x: 0.001044 ft^4
Iy-y: 0.003519 ft^4
J: 0.004563 ft^4
rx-x: 0.084304 ft
ry-y: 0.154797 ft
Sx-x (min): 0.005914 ft^3
Sy-y (min): 0.009830 ft^3
======================================================================
END OF REPORT
======================================================================
@stephensmitchell

Copy link
Copy Markdown
Author
# created with Alibre Script Genie by Stephen S. Mitchell, https://github.com/stephensmitchell
from __future__ import division
import sys
import math
from AlibreScript import *
ScriptName = "Batch Area Moments"
ScriptVersion = "1.0"
DIALOG_WIDTH = 400
CURVE_SEGMENTS = 360
VERTEX_TOLERANCE = 1e-6
CIRCUMFERENCE_TOLERANCE = 0.02
ARC_LENGTH_TOLERANCE = 0.05
CURVATURE_THRESHOLD = 1.05
def vec_subtract(a, b):
    return [a[0]-b[0], a[1]-b[1], a[2]-b[2]]
def vec_add(a, b):
    return [a[0]+b[0], a[1]+b[1], a[2]+b[2]]
def vec_scale(v, s):
    return [v[0]*s, v[1]*s, v[2]*s]
def vec_cross(a, b):
    return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]
def vec_dot(a, b):
    return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]
def vec_length(v):
    return math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
def vec_normalize(v):
    L = vec_length(v)
    if L < 1e-12:
        return [0.0, 0.0, 0.0]
    return [v[0]/L, v[1]/L, v[2]/L]
def vec_distance(a, b):
    return math.sqrt((b[0]-a[0])**2 + (b[1]-a[1])**2 + (b[2]-a[2])**2)
def vec_midpoint(a, b):
    return [(a[0]+b[0])/2.0, (a[1]+b[1])/2.0, (a[2]+b[2])/2.0]
def vec_lerp(a, b, t):
    return [a[0] + t*(b[0]-a[0]), a[1] + t*(b[1]-a[1]), a[2] + t*(b[2]-a[2])]
def safe_sqrt(val):
    if val < 0:
        return 0.0
    return math.sqrt(val)
def safe_div(num, denom, default=float('inf')):
    if abs(denom) < 1e-12:
        return default
    return num / denom
class CircleProperties:
    def __init__(self, radius, center_x=0.0, center_y=0.0):
        self.shape_type = "Circle"
        self.radius = abs(radius)
        self.n = 0
        r = self.radius
        A = math.pi * r * r
        I = math.pi * r**4 / 4.0
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = I
        self.Iy_centroid = I
        self.Ixy_centroid = 0.0
        self.J_centroid = 2.0 * I
        self.Ix_origin = I + A * center_y**2
        self.Iy_origin = I + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = I
        self.I_min = I
        self.theta_principal = 0.0
        self.theta_principal_deg = 0.0
        self.rx = r / 2.0
        self.ry = r / 2.0
        self.rp = r / math.sqrt(2.0)
        self.c_top = r
        self.c_bottom = r
        self.c_right = r
        self.c_left = r
        S = safe_div(I, r, 0.0)
        self.Sx_top = S
        self.Sx_bottom = S
        self.Sy_right = S
        self.Sy_left = S
        self.Sx_min = S
        self.Sy_min = S
class AnnulusProperties:
    def __init__(self, outer_radius, inner_radius, center_x=0.0, center_y=0.0):
        self.shape_type = "Annulus"
        R = abs(outer_radius)
        r = abs(inner_radius)
        if R < r:
            R, r = r, R
        self.outer_radius = R
        self.inner_radius = r
        self.n = 0
        A = math.pi * (R**2 - r**2)
        I = math.pi * (R**4 - r**4) / 4.0
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = I
        self.Iy_centroid = I
        self.Ixy_centroid = 0.0
        self.J_centroid = 2.0 * I
        self.Ix_origin = I + A * center_y**2
        self.Iy_origin = I + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = I
        self.I_min = I
        self.theta_principal = 0.0
        self.theta_principal_deg = 0.0
        self.rx = safe_sqrt(safe_div(I, A, 0.0))
        self.ry = safe_sqrt(safe_div(I, A, 0.0))
        self.rp = safe_sqrt(safe_div(2.0 * I, A, 0.0))
        self.c_top = R
        self.c_bottom = R
        self.c_right = R
        self.c_left = R
        S = safe_div(I, R, 0.0)
        self.Sx_top = S
        self.Sx_bottom = S
        self.Sy_right = S
        self.Sy_left = S
        self.Sx_min = S
        self.Sy_min = S
class RectangleProperties:
    def __init__(self, width, height, center_x=0.0, center_y=0.0, rotation=0.0):
        self.shape_type = "Rectangle"
        self.width = abs(width)
        self.height = abs(height)
        self.rotation = rotation
        self.n = 4
        b = self.width
        h = self.height
        A = b * h
        Ix_local = b * h**3 / 12.0
        Iy_local = b**3 * h / 12.0
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        cos2 = math.cos(2.0 * rotation)
        sin2 = math.sin(2.0 * rotation)
        I_avg = (Ix_local + Iy_local) / 2.0
        I_diff = (Ix_local - Iy_local) / 2.0
        self.Ix_centroid = I_avg + I_diff * cos2
        self.Iy_centroid = I_avg - I_diff * cos2
        self.Ixy_centroid = -I_diff * sin2
        self.J_centroid = Ix_local + Iy_local
        self.Ix_origin = self.Ix_centroid + A * center_y**2
        self.Iy_origin = self.Iy_centroid + A * center_x**2
        self.Ixy_origin = self.Ixy_centroid + A * center_x * center_y
        self.I_max = max(Ix_local, Iy_local)
        self.I_min = min(Ix_local, Iy_local)
        self.theta_principal = rotation if Ix_local >= Iy_local else rotation + math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        cos_r = abs(math.cos(rotation))
        sin_r = abs(math.sin(rotation))
        self.c_top = (h * cos_r + b * sin_r) / 2.0
        self.c_bottom = self.c_top
        self.c_right = (b * cos_r + h * sin_r) / 2.0
        self.c_left = self.c_right
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right
class SemicircleProperties:
    def __init__(self, radius, center_x=0.0, center_y=0.0, rotation=0.0):
        self.shape_type = "Semicircle"
        self.radius = abs(radius)
        self.rotation = rotation
        self.n = 0
        r = self.radius
        A = math.pi * r**2 / 2.0
        y_bar = 4.0 * r / (3.0 * math.pi)
        Ix_local = (math.pi/8.0 - 8.0/(9.0*math.pi)) * r**4
        Iy_local = math.pi * r**4 / 8.0
        self.area = A
        cos_r = math.cos(rotation)
        sin_r = math.sin(rotation)
        self.centroid_x = center_x + y_bar * sin_r
        self.centroid_y = center_y + y_bar * cos_r
        cos2 = math.cos(2.0 * rotation)
        sin2 = math.sin(2.0 * rotation)
        I_avg = (Ix_local + Iy_local) / 2.0
        I_diff = (Ix_local - Iy_local) / 2.0
        self.Ix_centroid = I_avg + I_diff * cos2
        self.Iy_centroid = I_avg - I_diff * cos2
        self.Ixy_centroid = -I_diff * sin2
        self.J_centroid = Ix_local + Iy_local
        self.Ix_origin = self.Ix_centroid + A * self.centroid_y**2
        self.Iy_origin = self.Iy_centroid + A * self.centroid_x**2
        self.Ixy_origin = self.Ixy_centroid + A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix_local, Iy_local)
        self.I_min = min(Ix_local, Iy_local)
        self.theta_principal = rotation
        self.theta_principal_deg = math.degrees(rotation)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = r - y_bar
        self.c_bottom = y_bar
        self.c_right = r
        self.c_left = r
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)
class PolygonProperties:
    def __init__(self, vertices_2d, shape_name="Polygon"):
        self.shape_type = shape_name
        self.vertices = vertices_2d
        self.n = len(vertices_2d)
        if self.n < 3:
            raise ValueError("Polygon must have at least 3 vertices")
        self._compute_all()
    def _compute_all(self):
        verts = self.vertices
        n = self.n
        signed_area = 0.0
        for i in range(n):
            j = (i + 1) % n
            signed_area += verts[i][0] * verts[j][1]
            signed_area -= verts[j][0] * verts[i][1]
        signed_area *= 0.5
        self.area = abs(signed_area)
        if self.area < 1e-12:
            raise ValueError("Polygon has zero or negligible area")
        cx, cy = 0.0, 0.0
        for i in range(n):
            j = (i + 1) % n
            cross = verts[i][0] * verts[j][1] - verts[j][0] * verts[i][1]
            cx += (verts[i][0] + verts[j][0]) * cross
            cy += (verts[i][1] + verts[j][1]) * cross
        factor = 1.0 / (6.0 * signed_area) if abs(signed_area) > 1e-12 else 0.0
        self.centroid_x = cx * factor
        self.centroid_y = cy * factor
        Ix, Iy, Ixy = 0.0, 0.0, 0.0
        for i in range(n):
            j = (i + 1) % n
            x0, y0 = verts[i]
            x1, y1 = verts[j]
            cross = x0 * y1 - x1 * y0
            Ix += (y0*y0 + y0*y1 + y1*y1) * cross
            Iy += (x0*x0 + x0*x1 + x1*x1) * cross
            Ixy += (x0*y1 + 2.0*x0*y0 + 2.0*x1*y1 + x1*y0) * cross
        self.Ix_origin = abs(Ix / 12.0)
        self.Iy_origin = abs(Iy / 12.0)
        self.Ixy_origin = Ixy / 24.0 if signed_area > 0 else -Ixy / 24.0
        A = self.area
        self.Ix_centroid = abs(self.Ix_origin - A * self.centroid_y**2)
        self.Iy_centroid = abs(self.Iy_origin - A * self.centroid_x**2)
        self.Ixy_centroid = self.Ixy_origin - A * self.centroid_x * self.centroid_y
        self.J_centroid = self.Ix_centroid + self.Iy_centroid
        Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
        I_avg = (Ix_c + Iy_c) / 2.0
        I_diff = (Ix_c - Iy_c) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy_c**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        x_rel = [v[0] - self.centroid_x for v in verts]
        y_rel = [v[1] - self.centroid_y for v in verts]
        self.c_top = max(y_rel) if y_rel else 0.0
        self.c_bottom = abs(min(y_rel)) if y_rel else 0.0
        self.c_right = max(x_rel) if x_rel else 0.0
        self.c_left = abs(min(x_rel)) if x_rel else 0.0
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)
class RectangularTubingProperties:
    def __init__(self, outer_width, outer_height, inner_width, inner_height, center_x=0.0, center_y=0.0):
        self.shape_type = "Rectangular Tubing"
        self.outer_width = abs(outer_width)
        self.outer_height = abs(outer_height)
        self.inner_width = abs(inner_width)
        self.inner_height = abs(inner_height)
        self.n = 8
        B, H = self.outer_width, self.outer_height
        b, h = self.inner_width, self.inner_height
        A_outer = B * H
        A_inner = b * h
        A = A_outer - A_inner
        Ix_outer = B * H**3 / 12.0
        Iy_outer = B**3 * H / 12.0
        Ix_inner = b * h**3 / 12.0
        Iy_inner = b**3 * h / 12.0
        Ix = Ix_outer - Ix_inner
        Iy = Iy_outer - Iy_inner
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = B / 2.0
        self.c_left = B / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class EllipseProperties:
    def __init__(self, semi_major, semi_minor, center_x=0.0, center_y=0.0, rotation=0.0):
        self.shape_type = "Ellipse"
        a = abs(semi_major)
        b = abs(semi_minor)
        if a < b:
            a, b = b, a
        self.semi_major = a
        self.semi_minor = b
        self.rotation = rotation
        self.n = 0
        A = math.pi * a * b
        Ix_local = math.pi * a * b**3 / 4.0
        Iy_local = math.pi * a**3 * b / 4.0
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        cos2 = math.cos(2.0 * rotation)
        sin2 = math.sin(2.0 * rotation)
        I_avg = (Ix_local + Iy_local) / 2.0
        I_diff = (Ix_local - Iy_local) / 2.0
        self.Ix_centroid = I_avg + I_diff * cos2
        self.Iy_centroid = I_avg - I_diff * cos2
        self.Ixy_centroid = -I_diff * sin2
        self.J_centroid = Ix_local + Iy_local
        self.Ix_origin = self.Ix_centroid + A * center_y**2
        self.Iy_origin = self.Iy_centroid + A * center_x**2
        self.Ixy_origin = self.Ixy_centroid + A * center_x * center_y
        self.I_max = max(Ix_local, Iy_local)
        self.I_min = min(Ix_local, Iy_local)
        self.theta_principal = rotation if Iy_local >= Ix_local else rotation + math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = b
        self.c_bottom = b
        self.c_right = a
        self.c_left = a
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class EllipticalTubeProperties:
    def __init__(self, outer_major, outer_minor, inner_major, inner_minor, center_x=0.0, center_y=0.0):
        self.shape_type = "Elliptical Tube"
        a1 = abs(outer_major)
        b1 = abs(outer_minor)
        a2 = abs(inner_major)
        b2 = abs(inner_minor)
        self.outer_major = a1
        self.outer_minor = b1
        self.inner_major = a2
        self.inner_minor = b2
        self.n = 0
        A_outer = math.pi * a1 * b1
        A_inner = math.pi * a2 * b2
        A = A_outer - A_inner
        Ix_outer = math.pi * a1 * b1**3 / 4.0
        Iy_outer = math.pi * a1**3 * b1 / 4.0
        Ix_inner = math.pi * a2 * b2**3 / 4.0
        Iy_inner = math.pi * a2**3 * b2 / 4.0
        Ix = Ix_outer - Ix_inner
        Iy = Iy_outer - Iy_inner
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Iy >= Ix else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = b1
        self.c_bottom = b1
        self.c_right = a1
        self.c_left = a1
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class QuarterCircleProperties:
    def __init__(self, radius, center_x=0.0, center_y=0.0, quadrant=1):
        self.shape_type = "Quarter Circle"
        self.radius = abs(radius)
        self.quadrant = quadrant
        self.n = 0
        r = self.radius
        A = math.pi * r**2 / 4.0
        d = 4.0 * r / (3.0 * math.pi)
        Ix_local = (math.pi/16.0 - 4.0/(9.0*math.pi)) * r**4
        Iy_local = Ix_local
        self.area = A
        if quadrant == 1:
            self.centroid_x = center_x + d
            self.centroid_y = center_y + d
        elif quadrant == 2:
            self.centroid_x = center_x - d
            self.centroid_y = center_y + d
        elif quadrant == 3:
            self.centroid_x = center_x - d
            self.centroid_y = center_y - d
        else:
            self.centroid_x = center_x + d
            self.centroid_y = center_y - d
        self.Ix_centroid = Ix_local
        self.Iy_centroid = Iy_local
        self.Ixy_centroid = (1.0/8.0 - 4.0/(9.0*math.pi)) * r**4
        if quadrant in (2, 4):
            self.Ixy_centroid = -self.Ixy_centroid
        self.J_centroid = Ix_local + Iy_local
        self.Ix_origin = self.Ix_centroid + A * self.centroid_y**2
        self.Iy_origin = self.Iy_centroid + A * self.centroid_x**2
        self.Ixy_origin = self.Ixy_centroid + A * self.centroid_x * self.centroid_y
        Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
        I_avg = (Ix_c + Iy_c) / 2.0
        I_diff = (Ix_c - Iy_c) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy_c**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = r - d if quadrant in (1, 2) else d
        self.c_bottom = d if quadrant in (1, 2) else r - d
        self.c_right = r - d if quadrant in (1, 4) else d
        self.c_left = d if quadrant in (1, 4) else r - d
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class SlotProperties:
    def __init__(self, length, width, center_x=0.0, center_y=0.0, rotation=0.0):
        self.shape_type = "Slot"
        L = abs(length)
        W = abs(width)
        if L < W:
            L, W = W, L
        self.length = L
        self.width = W
        self.rotation = rotation
        self.n = 0
        r = W / 2.0
        rect_len = L - W
        A_rect = rect_len * W
        A_circle = math.pi * r**2
        A = A_rect + A_circle
        Ix_rect = rect_len * W**3 / 12.0
        Iy_rect = rect_len**3 * W / 12.0
        Ix_circle = math.pi * r**4 / 4.0
        Iy_semi = math.pi * r**4 / 8.0
        d_semi = rect_len / 2.0
        Iy_semis = 2.0 * (Iy_semi + (A_circle/2.0) * d_semi**2)
        Ix_local = Ix_rect + Ix_circle
        Iy_local = Iy_rect + Iy_semis
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        cos2 = math.cos(2.0 * rotation)
        sin2 = math.sin(2.0 * rotation)
        I_avg = (Ix_local + Iy_local) / 2.0
        I_diff = (Ix_local - Iy_local) / 2.0
        self.Ix_centroid = I_avg + I_diff * cos2
        self.Iy_centroid = I_avg - I_diff * cos2
        self.Ixy_centroid = -I_diff * sin2
        self.J_centroid = Ix_local + Iy_local
        self.Ix_origin = self.Ix_centroid + A * center_y**2
        self.Iy_origin = self.Iy_centroid + A * center_x**2
        self.Ixy_origin = self.Ixy_centroid + A * center_x * center_y
        self.I_max = max(Ix_local, Iy_local)
        self.I_min = min(Ix_local, Iy_local)
        self.theta_principal = rotation if Iy_local >= Ix_local else rotation + math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, A, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = W / 2.0
        self.c_bottom = W / 2.0
        self.c_right = L / 2.0
        self.c_left = L / 2.0
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class ISectionProperties:
    def __init__(self, height, flange_width, web_thickness, flange_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "I-Section"
        H = abs(height)
        B = abs(flange_width)
        tw = abs(web_thickness)
        tf = abs(flange_thickness)
        self.height = H
        self.flange_width = B
        self.web_thickness = tw
        self.flange_thickness = tf
        self.n = 12
        A_flanges = 2.0 * B * tf
        A_web = (H - 2.0 * tf) * tw
        A = A_flanges + A_web
        Ix_flanges = 2.0 * (B * tf**3 / 12.0 + B * tf * ((H - tf) / 2.0)**2)
        Ix_web = tw * (H - 2.0 * tf)**3 / 12.0
        Ix = Ix_flanges + Ix_web
        Iy_flanges = 2.0 * tf * B**3 / 12.0
        Iy_web = (H - 2.0 * tf) * tw**3 / 12.0
        Iy = Iy_flanges + Iy_web
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = B / 2.0
        self.c_left = B / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class TSectionProperties:
    def __init__(self, height, flange_width, web_thickness, flange_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "T-Section"
        H = abs(height)
        B = abs(flange_width)
        tw = abs(web_thickness)
        tf = abs(flange_thickness)
        self.height = H
        self.flange_width = B
        self.web_thickness = tw
        self.flange_thickness = tf
        self.n = 8
        A_flange = B * tf
        A_web = (H - tf) * tw
        A = A_flange + A_web
        y_flange = H - tf / 2.0
        y_web = (H - tf) / 2.0
        y_bar = (A_flange * y_flange + A_web * y_web) / A
        d_flange = y_flange - y_bar
        d_web = y_web - y_bar
        Ix_flange = B * tf**3 / 12.0 + A_flange * d_flange**2
        Ix_web = tw * (H - tf)**3 / 12.0 + A_web * d_web**2
        Ix = Ix_flange + Ix_web
        Iy_flange = tf * B**3 / 12.0
        Iy_web = (H - tf) * tw**3 / 12.0
        Iy = Iy_flange + Iy_web
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        self.c_right = B / 2.0
        self.c_left = B / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = self.Sy_right

class ChannelSectionProperties:
    def __init__(self, height, flange_width, web_thickness, flange_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Channel Section"
        H = abs(height)
        B = abs(flange_width)
        tw = abs(web_thickness)
        tf = abs(flange_thickness)
        self.height = H
        self.flange_width = B
        self.web_thickness = tw
        self.flange_thickness = tf
        self.n = 8
        A_flanges = 2.0 * B * tf
        A_web = (H - 2.0 * tf) * tw
        A = A_flanges + A_web
        x_flanges = B / 2.0
        x_web = tw / 2.0
        x_bar = (A_flanges * x_flanges + A_web * x_web) / A
        Ix_flanges = 2.0 * (B * tf**3 / 12.0 + B * tf * ((H - tf) / 2.0)**2)
        Ix_web = tw * (H - 2.0 * tf)**3 / 12.0
        Ix = Ix_flanges + Ix_web
        d_flange = x_flanges - x_bar
        d_web = x_web - x_bar
        Iy_flanges = 2.0 * (tf * B**3 / 12.0 + B * tf * d_flange**2)
        Iy_web = (H - 2.0 * tf) * tw**3 / 12.0 + A_web * d_web**2
        Iy = Iy_flanges + Iy_web
        self.area = A
        self.centroid_x = center_x + (x_bar - B / 2.0)
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = B - x_bar
        self.c_left = x_bar
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = safe_div(Iy, self.c_left)
        self.Sx_min = self.Sx_top
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class AngleSectionProperties:
    def __init__(self, leg_height, leg_width, thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Angle Section"
        H = abs(leg_height)
        B = abs(leg_width)
        t = abs(thickness)
        self.leg_height = H
        self.leg_width = B
        self.thickness = t
        self.n = 6
        A_vert = H * t
        A_horiz = (B - t) * t
        A = A_vert + A_horiz
        x_vert = t / 2.0
        x_horiz = t + (B - t) / 2.0
        y_vert = H / 2.0
        y_horiz = t / 2.0
        x_bar = (A_vert * x_vert + A_horiz * x_horiz) / A
        y_bar = (A_vert * y_vert + A_horiz * y_horiz) / A
        dx_vert = x_vert - x_bar
        dy_vert = y_vert - y_bar
        dx_horiz = x_horiz - x_bar
        dy_horiz = y_horiz - y_bar
        Ix_vert = t * H**3 / 12.0 + A_vert * dy_vert**2
        Ix_horiz = (B - t) * t**3 / 12.0 + A_horiz * dy_horiz**2
        Ix = Ix_vert + Ix_horiz
        Iy_vert = H * t**3 / 12.0 + A_vert * dx_vert**2
        Iy_horiz = t * (B - t)**3 / 12.0 + A_horiz * dx_horiz**2
        Iy = Iy_vert + Iy_horiz
        Ixy_vert = A_vert * dx_vert * dy_vert
        Ixy_horiz = A_horiz * dx_horiz * dy_horiz
        Ixy = Ixy_vert + Ixy_horiz
        self.area = A
        self.centroid_x = center_x + (x_bar - B / 2.0)
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = Ixy
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = Ixy + A * self.centroid_x * self.centroid_y
        I_avg = (Ix + Iy) / 2.0
        I_diff = (Ix - Iy) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy, Ix - Iy)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        self.c_right = B - x_bar
        self.c_left = x_bar
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = safe_div(Iy, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class ZSectionProperties:
    def __init__(self, height, flange_width, web_thickness, flange_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Z-Section"
        H = abs(height)
        B = abs(flange_width)
        tw = abs(web_thickness)
        tf = abs(flange_thickness)
        self.height = H
        self.flange_width = B
        self.web_thickness = tw
        self.flange_thickness = tf
        self.n = 8
        A_top = B * tf
        A_bot = B * tf
        A_web = (H - 2.0 * tf) * tw
        A = A_top + A_bot + A_web
        Ix_top = B * tf**3 / 12.0 + A_top * ((H - tf) / 2.0)**2
        Ix_bot = B * tf**3 / 12.0 + A_bot * ((H - tf) / 2.0)**2
        Ix_web = tw * (H - 2.0 * tf)**3 / 12.0
        Ix = Ix_top + Ix_bot + Ix_web
        x_top = (B - tw) / 2.0
        x_bot = -(B - tw) / 2.0
        Iy_top = tf * B**3 / 12.0 + A_top * x_top**2
        Iy_bot = tf * B**3 / 12.0 + A_bot * x_bot**2
        Iy_web = (H - 2.0 * tf) * tw**3 / 12.0
        Iy = Iy_top + Iy_bot + Iy_web
        Ixy_top = A_top * x_top * ((H - tf) / 2.0)
        Ixy_bot = A_bot * x_bot * (-(H - tf) / 2.0)
        Ixy = Ixy_top + Ixy_bot
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = Ixy
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = Ixy + A * center_x * center_y
        I_avg = (Ix + Iy) / 2.0
        I_diff = (Ix - Iy) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy, Ix - Iy)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = B / 2.0 + (B - tw) / 2.0
        self.c_left = B / 2.0 + (B - tw) / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class HatSectionProperties:
    def __init__(self, height, top_width, bottom_width, thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Hat Section"
        H = abs(height)
        Bt = abs(top_width)
        Bb = abs(bottom_width)
        t = abs(thickness)
        self.height = H
        self.top_width = Bt
        self.bottom_width = Bb
        self.thickness = t
        self.n = 8
        A_top = Bt * t
        A_legs = 2.0 * (H - t) * t
        A_flanges = 2.0 * ((Bb - Bt) / 2.0) * t
        A = A_top + A_legs + A_flanges
        y_top = H - t / 2.0
        y_legs = (H - t) / 2.0 + t
        y_flanges = t / 2.0
        y_bar = (A_top * y_top + A_legs * y_legs + A_flanges * y_flanges) / A
        d_top = y_top - y_bar
        d_legs = y_legs - y_bar
        d_flanges = y_flanges - y_bar
        Ix_top = Bt * t**3 / 12.0 + A_top * d_top**2
        Ix_legs = 2.0 * (t * (H - t)**3 / 12.0 + (H - t) * t * d_legs**2)
        Ix_flanges = 2.0 * (((Bb - Bt) / 2.0) * t**3 / 12.0 + ((Bb - Bt) / 2.0) * t * d_flanges**2)
        Ix = Ix_top + Ix_legs + Ix_flanges
        Iy_top = t * Bt**3 / 12.0
        x_leg = Bt / 2.0 + t / 2.0
        Iy_legs = 2.0 * ((H - t) * t**3 / 12.0 + (H - t) * t * x_leg**2)
        x_flange = Bt / 2.0 + t + (Bb - Bt - 2.0 * t) / 4.0
        Iy_flanges = 2.0 * (t * ((Bb - Bt) / 2.0)**3 / 12.0 + ((Bb - Bt) / 2.0) * t * x_flange**2)
        Iy = Iy_top + Iy_legs + Iy_flanges
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        self.c_right = Bb / 2.0
        self.c_left = Bb / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = self.Sy_right

class SigmaSectionProperties:
    def __init__(self, height, flange_width, lip_height, thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Sigma Section"
        H = abs(height)
        B = abs(flange_width)
        L = abs(lip_height)
        t = abs(thickness)
        self.height = H
        self.flange_width = B
        self.lip_height = L
        self.thickness = t
        self.n = 12
        A_lips = 2.0 * L * t
        A_flanges = 2.0 * B * t
        A_web = (H - 2.0 * (L + t)) * t
        A = A_lips + A_flanges + A_web
        y_lip_top = H - L / 2.0
        y_lip_bot = L / 2.0
        y_flange_top = H - L - t / 2.0
        y_flange_bot = L + t / 2.0
        y_web = H / 2.0
        Ix = 0.0
        Ix += 2.0 * (t * L**3 / 12.0)
        Ix += L * t * (H - L / 2.0 - H / 2.0)**2
        Ix += L * t * (L / 2.0 - H / 2.0)**2
        Ix += 2.0 * (B * t**3 / 12.0)
        Ix += B * t * (H - L - t / 2.0 - H / 2.0)**2
        Ix += B * t * (L + t / 2.0 - H / 2.0)**2
        Ix += t * (H - 2.0 * (L + t))**3 / 12.0
        x_lips = B + t / 2.0
        x_flanges = B / 2.0
        x_web = t / 2.0
        x_bar = (A_lips * x_lips + A_flanges * x_flanges + A_web * x_web) / A
        d_lips = x_lips - x_bar
        d_flanges = x_flanges - x_bar
        d_web = x_web - x_bar
        Iy = 0.0
        Iy += 2.0 * (L * t**3 / 12.0 + L * t * d_lips**2)
        Iy += 2.0 * (t * B**3 / 12.0 + B * t * d_flanges**2)
        Iy += (H - 2.0 * (L + t)) * t**3 / 12.0 + A_web * d_web**2
        self.area = A
        self.centroid_x = center_x + (x_bar - (B + t) / 2.0)
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = B + t - x_bar
        self.c_left = x_bar
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = safe_div(Iy, self.c_left)
        self.Sx_min = self.Sx_top
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class OmegaSectionProperties:
    def __init__(self, height, top_width, leg_width, thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Omega Section"
        H = abs(height)
        Bt = abs(top_width)
        Bl = abs(leg_width)
        t = abs(thickness)
        self.height = H
        self.top_width = Bt
        self.leg_width = Bl
        self.thickness = t
        self.n = 10
        A_top = Bt * t
        A_legs = 2.0 * (H - t) * t
        A_feet = 2.0 * Bl * t
        A = A_top + A_legs + A_feet
        y_top = H - t / 2.0
        y_legs = t + (H - t) / 2.0
        y_feet = t / 2.0
        y_bar = (A_top * y_top + A_legs * y_legs + A_feet * y_feet) / A
        d_top = y_top - y_bar
        d_legs = y_legs - y_bar
        d_feet = y_feet - y_bar
        Ix_top = Bt * t**3 / 12.0 + A_top * d_top**2
        Ix_legs = 2.0 * (t * (H - t)**3 / 12.0 + (H - t) * t * d_legs**2)
        Ix_feet = 2.0 * (Bl * t**3 / 12.0 + Bl * t * d_feet**2)
        Ix = Ix_top + Ix_legs + Ix_feet
        Iy_top = t * Bt**3 / 12.0
        x_leg = Bt / 2.0 + t / 2.0
        Iy_legs = 2.0 * ((H - t) * t**3 / 12.0 + (H - t) * t * x_leg**2)
        x_foot = Bt / 2.0 + t + Bl / 2.0
        Iy_feet = 2.0 * (t * Bl**3 / 12.0 + Bl * t * x_foot**2)
        Iy = Iy_top + Iy_legs + Iy_feet
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        total_width = Bt + 2.0 * t + 2.0 * Bl
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        self.c_right = total_width / 2.0
        self.c_left = total_width / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = self.Sy_right

class DoubleAngleProperties:
    def __init__(self, leg_height, leg_width, thickness, gap, center_x=0.0, center_y=0.0):
        self.shape_type = "Double Angle"
        H = abs(leg_height)
        B = abs(leg_width)
        t = abs(thickness)
        g = abs(gap)
        self.leg_height = H
        self.leg_width = B
        self.thickness = t
        self.gap = g
        self.n = 12
        single = AngleSectionProperties(H, B, t, 0.0, 0.0)
        A = 2.0 * single.area
        Ix = 2.0 * single.Ix_centroid
        d_y = g / 2.0 + (B - single.centroid_x)
        Iy = 2.0 * (single.Iy_centroid + single.area * d_y**2)
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y + single.centroid_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = single.c_top
        self.c_bottom = single.c_bottom
        self.c_right = g / 2.0 + B
        self.c_left = self.c_right
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = self.Sy_right

class DoubleChannelProperties:
    def __init__(self, height, flange_width, web_thickness, flange_thickness, gap, center_x=0.0, center_y=0.0):
        self.shape_type = "Double Channel"
        H = abs(height)
        B = abs(flange_width)
        tw = abs(web_thickness)
        tf = abs(flange_thickness)
        g = abs(gap)
        self.height = H
        self.flange_width = B
        self.web_thickness = tw
        self.flange_thickness = tf
        self.gap = g
        self.n = 16
        single = ChannelSectionProperties(H, B, tw, tf, 0.0, 0.0)
        A = 2.0 * single.area
        Ix = 2.0 * single.Ix_centroid
        d_x = g / 2.0 + (B - abs(single.centroid_x))
        Iy = 2.0 * (single.Iy_centroid + single.area * d_x**2)
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = g / 2.0 + B
        self.c_left = self.c_right
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right

class BuiltUpISectionProperties:
    def __init__(self, height, top_flange_width, top_flange_thickness, bottom_flange_width, bottom_flange_thickness, web_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Built-Up I-Section"
        H = abs(height)
        Bt = abs(top_flange_width)
        tt = abs(top_flange_thickness)
        Bb = abs(bottom_flange_width)
        tb = abs(bottom_flange_thickness)
        tw = abs(web_thickness)
        self.height = H
        self.top_flange_width = Bt
        self.top_flange_thickness = tt
        self.bottom_flange_width = Bb
        self.bottom_flange_thickness = tb
        self.web_thickness = tw
        self.n = 12
        A_top = Bt * tt
        A_bot = Bb * tb
        A_web = (H - tt - tb) * tw
        A = A_top + A_bot + A_web
        y_top = H - tt / 2.0
        y_bot = tb / 2.0
        y_web = tb + (H - tt - tb) / 2.0
        y_bar = (A_top * y_top + A_bot * y_bot + A_web * y_web) / A
        d_top = y_top - y_bar
        d_bot = y_bot - y_bar
        d_web = y_web - y_bar
        Ix_top = Bt * tt**3 / 12.0 + A_top * d_top**2
        Ix_bot = Bb * tb**3 / 12.0 + A_bot * d_bot**2
        Ix_web = tw * (H - tt - tb)**3 / 12.0 + A_web * d_web**2
        Ix = Ix_top + Ix_bot + Ix_web
        Iy_top = tt * Bt**3 / 12.0
        Iy_bot = tb * Bb**3 / 12.0
        Iy_web = (H - tt - tb) * tw**3 / 12.0
        Iy = Iy_top + Iy_bot + Iy_web
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        max_flange = max(Bt, Bb)
        self.c_right = max_flange / 2.0
        self.c_left = max_flange / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = self.Sy_right

class BuiltUpBoxSectionProperties:
    def __init__(self, height, width, top_thickness, bottom_thickness, left_thickness, right_thickness, center_x=0.0, center_y=0.0):
        self.shape_type = "Built-Up Box Section"
        H = abs(height)
        W = abs(width)
        tt = abs(top_thickness)
        tb = abs(bottom_thickness)
        tl = abs(left_thickness)
        tr = abs(right_thickness)
        self.height = H
        self.width = W
        self.top_thickness = tt
        self.bottom_thickness = tb
        self.left_thickness = tl
        self.right_thickness = tr
        self.n = 8
        A_outer = H * W
        inner_h = H - tt - tb
        inner_w = W - tl - tr
        A_inner = inner_h * inner_w
        A = A_outer - A_inner
        y_outer = H / 2.0
        y_inner = tb + inner_h / 2.0
        x_outer = W / 2.0
        x_inner = tl + inner_w / 2.0
        y_bar = (A_outer * y_outer - A_inner * y_inner) / A
        x_bar = (A_outer * x_outer - A_inner * x_inner) / A
        Ix_outer = W * H**3 / 12.0
        Ix_inner = inner_w * inner_h**3 / 12.0
        d_outer = y_outer - y_bar
        d_inner = y_inner - y_bar
        Ix = Ix_outer + A_outer * d_outer**2 - (Ix_inner + A_inner * d_inner**2)
        Iy_outer = H * W**3 / 12.0
        Iy_inner = inner_h * inner_w**3 / 12.0
        dx_outer = x_outer - x_bar
        dx_inner = x_inner - x_bar
        Iy = Iy_outer + A_outer * dx_outer**2 - (Iy_inner + A_inner * dx_inner**2)
        self.area = A
        self.centroid_x = center_x + (x_bar - W / 2.0)
        self.centroid_y = center_y + (y_bar - H / 2.0)
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * self.centroid_y**2
        self.Iy_origin = Iy + A * self.centroid_x**2
        self.Ixy_origin = A * self.centroid_x * self.centroid_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H - y_bar
        self.c_bottom = y_bar
        self.c_right = W - x_bar
        self.c_left = x_bar
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = safe_div(Iy, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class PolygonWithPolygonalHoleProperties:
    def __init__(self, outer_vertices_2d, inner_vertices_2d, shape_name="Polygonal Tube"):
        self.shape_type = shape_name
        outer = PolygonProperties(outer_vertices_2d, "outer")
        inner = PolygonProperties(inner_vertices_2d, "inner")
        dx = inner.centroid_x - outer.centroid_x
        dy = inner.centroid_y - outer.centroid_y
        inner_Ix_at_centroid = inner.Ix_centroid + inner.area * dy**2
        inner_Iy_at_centroid = inner.Iy_centroid + inner.area * dx**2
        inner_Ixy_at_centroid = inner.Ixy_centroid + inner.area * dx * dy
        self.n = outer.n + inner.n
        self.area = outer.area - inner.area
        if self.area < 1e-12:
            raise ValueError("Inner area exceeds outer area")
        self.centroid_x = outer.centroid_x
        self.centroid_y = outer.centroid_y
        self.Ix_centroid = outer.Ix_centroid - inner_Ix_at_centroid
        self.Iy_centroid = outer.Iy_centroid - inner_Iy_at_centroid
        self.Ixy_centroid = outer.Ixy_centroid - inner_Ixy_at_centroid
        self.J_centroid = self.Ix_centroid + self.Iy_centroid
        self.Ix_origin = outer.Ix_origin - inner.Ix_origin
        self.Iy_origin = outer.Iy_origin - inner.Iy_origin
        self.Ixy_origin = outer.Ixy_origin - inner.Ixy_origin
        Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
        I_avg = (Ix_c + Iy_c) / 2.0
        I_diff = (Ix_c - Iy_c) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy_c**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, self.area, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, self.area, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, self.area, 0.0))
        self.c_top = outer.c_top
        self.c_bottom = outer.c_bottom
        self.c_right = outer.c_right
        self.c_left = outer.c_left
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)
        self.outer_vertices = outer_vertices_2d
        self.inner_vertices = inner_vertices_2d

class PerforatedPlateProperties:
    def __init__(self, width, height, hole_diameter, hole_spacing_x, hole_spacing_y, center_x=0.0, center_y=0.0):
        self.shape_type = "Perforated Plate"
        W = abs(width)
        H = abs(height)
        d = abs(hole_diameter)
        sx = abs(hole_spacing_x)
        sy = abs(hole_spacing_y)
        self.width = W
        self.height = H
        self.hole_diameter = d
        self.hole_spacing_x = sx
        self.hole_spacing_y = sy
        self.n = 4
        A_plate = W * H
        nx = max(1, int(W / sx)) if sx > 0 else 1
        ny = max(1, int(H / sy)) if sy > 0 else 1
        n_holes = nx * ny
        A_holes = n_holes * math.pi * (d / 2.0)**2
        A = A_plate - A_holes
        Ix_plate = W * H**3 / 12.0
        Iy_plate = H * W**3 / 12.0
        Ix_hole = math.pi * (d / 2.0)**4 / 4.0
        Iy_hole = Ix_hole
        Ix_holes = n_holes * Ix_hole
        Iy_holes = n_holes * Iy_hole
        Ix = Ix_plate - Ix_holes
        Iy = Iy_plate - Iy_holes
        self.area = A
        self.centroid_x = center_x
        self.centroid_y = center_y
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = 0.0
        self.J_centroid = Ix + Iy
        self.Ix_origin = Ix + A * center_y**2
        self.Iy_origin = Iy + A * center_x**2
        self.Ixy_origin = A * center_x * center_y
        self.I_max = max(Ix, Iy)
        self.I_min = min(Ix, Iy)
        self.theta_principal = 0.0 if Ix >= Iy else math.pi/2.0
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, A, 0.0))
        self.c_top = H / 2.0
        self.c_bottom = H / 2.0
        self.c_right = W / 2.0
        self.c_left = W / 2.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = self.Sx_top
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = self.Sy_right
        self.Sx_min = self.Sx_top
        self.Sy_min = self.Sy_right
        self.num_holes = n_holes

class ShipSectionProperties:
    def __init__(self, outer_boundary_2d, inner_boundaries_2d=None, baseline_y=None, shape_name="Ship Section"):
        self.shape_type = shape_name
        if inner_boundaries_2d is None:
            inner_boundaries_2d = []
        outer = PolygonProperties(outer_boundary_2d, "outer")
        total_hole_area = 0.0
        total_hole_Ix = 0.0
        total_hole_Iy = 0.0
        total_hole_Ixy = 0.0
        hole_Ax_sum = 0.0
        hole_Ay_sum = 0.0
        self.num_openings = len(inner_boundaries_2d)
        self.inner_boundaries = inner_boundaries_2d
        for inner_verts in inner_boundaries_2d:
            if len(inner_verts) >= 3:
                try:
                    inner = PolygonProperties(inner_verts, "inner")
                    total_hole_area += inner.area
                    dx = inner.centroid_x - outer.centroid_x
                    dy = inner.centroid_y - outer.centroid_y
                    total_hole_Ix += inner.Ix_centroid + inner.area * dy**2
                    total_hole_Iy += inner.Iy_centroid + inner.area * dx**2
                    total_hole_Ixy += inner.Ixy_centroid + inner.area * dx * dy
                    hole_Ax_sum += inner.area * inner.centroid_x
                    hole_Ay_sum += inner.area * inner.centroid_y
                except Exception:
                    pass
        net_area = outer.area - total_hole_area
        if net_area < 1e-12:
            raise ValueError("Total hole area exceeds outer boundary area")
        net_Ax = outer.area * outer.centroid_x - hole_Ax_sum
        net_Ay = outer.area * outer.centroid_y - hole_Ay_sum
        net_cx = net_Ax / net_area
        net_cy = net_Ay / net_area
        dx_outer = outer.centroid_x - net_cx
        dy_outer = outer.centroid_y - net_cy
        Ix_outer_shifted = outer.Ix_centroid + outer.area * dy_outer**2
        Iy_outer_shifted = outer.Iy_centroid + outer.area * dx_outer**2
        Ixy_outer_shifted = outer.Ixy_centroid + outer.area * dx_outer * dy_outer
        total_hole_Ix_shifted = 0.0
        total_hole_Iy_shifted = 0.0
        total_hole_Ixy_shifted = 0.0
        for inner_verts in inner_boundaries_2d:
            if len(inner_verts) >= 3:
                try:
                    inner = PolygonProperties(inner_verts, "inner")
                    dx = inner.centroid_x - net_cx
                    dy = inner.centroid_y - net_cy
                    total_hole_Ix_shifted += inner.Ix_centroid + inner.area * dy**2
                    total_hole_Iy_shifted += inner.Iy_centroid + inner.area * dx**2
                    total_hole_Ixy_shifted += inner.Ixy_centroid + inner.area * dx * dy
                except Exception:
                    pass
        self.area = net_area
        self.centroid_x = net_cx
        self.centroid_y = net_cy
        self.Ix_centroid = Ix_outer_shifted - total_hole_Ix_shifted
        self.Iy_centroid = Iy_outer_shifted - total_hole_Iy_shifted
        self.Ixy_centroid = Ixy_outer_shifted - total_hole_Ixy_shifted
        self.J_centroid = self.Ix_centroid + self.Iy_centroid
        self.Ix_origin = outer.Ix_origin
        for inner_verts in inner_boundaries_2d:
            if len(inner_verts) >= 3:
                try:
                    inner = PolygonProperties(inner_verts, "inner")
                    self.Ix_origin -= inner.Ix_origin
                except Exception:
                    pass
        self.Iy_origin = outer.Iy_origin
        for inner_verts in inner_boundaries_2d:
            if len(inner_verts) >= 3:
                try:
                    inner = PolygonProperties(inner_verts, "inner")
                    self.Iy_origin -= inner.Iy_origin
                except Exception:
                    pass
        self.Ixy_origin = outer.Ixy_origin
        for inner_verts in inner_boundaries_2d:
            if len(inner_verts) >= 3:
                try:
                    inner = PolygonProperties(inner_verts, "inner")
                    self.Ixy_origin -= inner.Ixy_origin
                except Exception:
                    pass
        Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
        I_avg = (Ix_c + Iy_c) / 2.0
        I_diff = (Ix_c - Iy_c) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy_c**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, self.area, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, self.area, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, self.area, 0.0))
        all_y = [v[1] for v in outer_boundary_2d]
        all_x = [v[0] for v in outer_boundary_2d]
        y_max = max(all_y)
        y_min = min(all_y)
        x_max = max(all_x)
        x_min = min(all_x)
        self.height = y_max - y_min
        self.width = x_max - x_min
        self.c_top = y_max - net_cy
        self.c_bottom = net_cy - y_min
        self.c_right = x_max - net_cx
        self.c_left = net_cx - x_min
        if baseline_y is None:
            baseline_y = y_min
        self.baseline_y = baseline_y
        self.neutral_axis_height = net_cy - baseline_y
        self.c_deck = y_max - net_cy
        self.c_keel = net_cy - baseline_y
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_deck = safe_div(self.Ix_centroid, self.c_deck)
        self.Sx_keel = safe_div(self.Ix_centroid, self.c_keel)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)
        self.n = outer.n + sum(len(b) for b in inner_boundaries_2d if len(b) >= 3)

class CompositeRegionProperties:
    def __init__(self, regions, shape_name="Composite Region"):
        self.shape_type = shape_name
        self.regions = regions
        self.n = sum(r.n for r in regions if hasattr(r, 'n'))
        total_A = sum(r.area for r in regions)
        if total_A < 1e-12:
            raise ValueError("Composite has zero area")
        cx = sum(r.area * r.centroid_x for r in regions) / total_A
        cy = sum(r.area * r.centroid_y for r in regions) / total_A
        self.area = total_A
        self.centroid_x = cx
        self.centroid_y = cy
        Ix = 0.0
        Iy = 0.0
        Ixy = 0.0
        for r in regions:
            dx = r.centroid_x - cx
            dy = r.centroid_y - cy
            Ix += r.Ix_centroid + r.area * dy**2
            Iy += r.Iy_centroid + r.area * dx**2
            Ixy += r.Ixy_centroid + r.area * dx * dy
        self.Ix_centroid = Ix
        self.Iy_centroid = Iy
        self.Ixy_centroid = Ixy
        self.J_centroid = Ix + Iy
        self.Ix_origin = sum(r.Ix_origin for r in regions)
        self.Iy_origin = sum(r.Iy_origin for r in regions)
        self.Ixy_origin = sum(r.Ixy_origin for r in regions)
        I_avg = (Ix + Iy) / 2.0
        I_diff = (Ix - Iy) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy, Ix - Iy)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(Ix, total_A, 0.0))
        self.ry = safe_sqrt(safe_div(Iy, total_A, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, total_A, 0.0))
        all_c_top = [r.c_top + (r.centroid_y - cy) for r in regions if hasattr(r, 'c_top')]
        all_c_bot = [r.c_bottom + (cy - r.centroid_y) for r in regions if hasattr(r, 'c_bottom')]
        all_c_right = [r.c_right + (r.centroid_x - cx) for r in regions if hasattr(r, 'c_right')]
        all_c_left = [r.c_left + (cx - r.centroid_x) for r in regions if hasattr(r, 'c_left')]
        self.c_top = max(all_c_top) if all_c_top else 0.0
        self.c_bottom = max(all_c_bot) if all_c_bot else 0.0
        self.c_right = max(all_c_right) if all_c_right else 0.0
        self.c_left = max(all_c_left) if all_c_left else 0.0
        self.Sx_top = safe_div(Ix, self.c_top)
        self.Sx_bottom = safe_div(Ix, self.c_bottom)
        self.Sy_right = safe_div(Iy, self.c_right)
        self.Sy_left = safe_div(Iy, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)

class PolygonWithCircularHoleProperties:
    def __init__(self, outer_vertices_2d, hole_radius, hole_center_x=None, hole_center_y=None, shape_name="Polygon with Hole"):
        self.shape_type = shape_name
        outer = PolygonProperties(outer_vertices_2d, "outer")
        if hole_center_x is None:
            hole_center_x = outer.centroid_x
        if hole_center_y is None:
            hole_center_y = outer.centroid_y
        r = abs(hole_radius)
        hole_area = math.pi * r * r
        hole_Ic = math.pi * r**4 / 4.0
        dx = hole_center_x - outer.centroid_x
        dy = hole_center_y - outer.centroid_y
        hole_Ix_at_centroid = hole_Ic + hole_area * dy**2
        hole_Iy_at_centroid = hole_Ic + hole_area * dx**2
        hole_Ixy_at_centroid = hole_area * dx * dy
        self.n = outer.n
        self.area = outer.area - hole_area
        if self.area < 1e-12:
            raise ValueError("Hole area exceeds outer polygon area")
        self.centroid_x = outer.centroid_x
        self.centroid_y = outer.centroid_y
        self.Ix_centroid = outer.Ix_centroid - hole_Ix_at_centroid
        self.Iy_centroid = outer.Iy_centroid - hole_Iy_at_centroid
        self.Ixy_centroid = outer.Ixy_centroid - hole_Ixy_at_centroid
        self.J_centroid = self.Ix_centroid + self.Iy_centroid
        self.Ix_origin = outer.Ix_origin - (hole_Ic + hole_area * hole_center_y**2)
        self.Iy_origin = outer.Iy_origin - (hole_Ic + hole_area * hole_center_x**2)
        self.Ixy_origin = outer.Ixy_origin - hole_area * hole_center_x * hole_center_y
        Ix_c, Iy_c, Ixy_c = self.Ix_centroid, self.Iy_centroid, self.Ixy_centroid
        I_avg = (Ix_c + Iy_c) / 2.0
        I_diff = (Ix_c - Iy_c) / 2.0
        R = safe_sqrt(I_diff**2 + Ixy_c**2)
        self.I_max = I_avg + R
        self.I_min = max(0.0, I_avg - R)
        if abs(Ixy_c) < 1e-12 and abs(I_diff) < 1e-12:
            self.theta_principal = 0.0
        else:
            self.theta_principal = 0.5 * math.atan2(-2.0 * Ixy_c, Ix_c - Iy_c)
        self.theta_principal_deg = math.degrees(self.theta_principal)
        self.rx = safe_sqrt(safe_div(self.Ix_centroid, self.area, 0.0))
        self.ry = safe_sqrt(safe_div(self.Iy_centroid, self.area, 0.0))
        self.rp = safe_sqrt(safe_div(self.J_centroid, self.area, 0.0))
        self.c_top = outer.c_top
        self.c_bottom = outer.c_bottom
        self.c_right = outer.c_right
        self.c_left = outer.c_left
        self.Sx_top = safe_div(self.Ix_centroid, self.c_top)
        self.Sx_bottom = safe_div(self.Ix_centroid, self.c_bottom)
        self.Sy_right = safe_div(self.Iy_centroid, self.c_right)
        self.Sy_left = safe_div(self.Iy_centroid, self.c_left)
        self.Sx_min = min(self.Sx_top, self.Sx_bottom)
        self.Sy_min = min(self.Sy_right, self.Sy_left)
        self.outer_area = outer.area
        self.hole_radius = r
        self.hole_area = hole_area
        x_coords = [v[0] for v in outer_vertices_2d]
        y_coords = [v[1] for v in outer_vertices_2d]
        self.outer_width = max(x_coords) - min(x_coords)
        self.outer_height = max(y_coords) - min(y_coords)
class EdgeInfo:
    def __init__(self, edge, index):
        self.edge = edge
        self.index = index
        self.diameter = None
        self.length = None
        self.vertices = []
        self.vertex_count = 0
        self.is_circular = False
        self.is_closed = False
        self.is_full_circle = False
        self.is_curved = False
        self.chord_length = 0.0
        self._analyze()
    def _analyze(self):
        try:
            d = self.edge.Diameter
            if d is not None and d > 0:
                self.diameter = d
                self.is_circular = True
        except Exception:
            pass
        try:
            self.length = self.edge.Length
        except Exception:
            pass
        try:
            verts = self.edge.GetVertices()
            if verts:
                self.vertex_count = len(verts)
                for v in verts:
                    self.vertices.append([v.X, v.Y, v.Z])
        except Exception:
            pass
        self.is_closed = (self.vertex_count == 0)
        if self.is_circular and self.diameter:
            if self.is_closed:
                self.is_full_circle = True
            elif not self.length or self.length < 1e-9:
                self.is_full_circle = True
            elif self.length > 0:
                expected_circumference = math.pi * self.diameter
                if abs(self.length - expected_circumference) < expected_circumference * CIRCUMFERENCE_TOLERANCE:
                    self.is_full_circle = True
        if self.vertex_count >= 2:
            self.chord_length = vec_distance(self.vertices[0], self.vertices[1])
            if self.length and self.chord_length > 1e-12:
                ratio = self.length / self.chord_length
                if ratio > CURVATURE_THRESHOLD:
                    self.is_curved = True
def is_complete_circle(edge_info):
    return edge_info.is_closed or edge_info.is_full_circle

def find_edge_loops(edge_infos):
    if not edge_infos:
        return []
    remaining = list(edge_infos)
    loops = []
    while remaining:
        loop = [remaining.pop(0)]
        changed = True
        while changed:
            changed = False
            for i, edge in enumerate(remaining):
                if edges_share_vertex(loop[-1], edge):
                    loop.append(remaining.pop(i))
                    changed = True
                    break
                elif edges_share_vertex(loop[0], edge):
                    loop.insert(0, remaining.pop(i))
                    changed = True
                    break
        loops.append(loop)
    return loops

def get_rectangle_dimensions(loop_edges):
    if len(loop_edges) != 4:
        return None
    lengths = []
    for e in loop_edges:
        if e.length is None or e.length < 1e-9:
            return None
        lengths.append(e.length)
    lengths_sorted = sorted(lengths)
    if abs(lengths_sorted[0] - lengths_sorted[1]) > 0.01:
        return None
    if abs(lengths_sorted[2] - lengths_sorted[3]) > 0.01:
        return None
    width = (lengths_sorted[2] + lengths_sorted[3]) / 2.0
    height = (lengths_sorted[0] + lengths_sorted[1]) / 2.0
    return (width, height)
def collect_unique_vertices(edges, tolerance=VERTEX_TOLERANCE):
    all_verts = []
    for e in edges:
        for v in e.vertices:
            is_dup = False
            for existing in all_verts:
                if vec_distance(v, existing) < tolerance:
                    is_dup = True
                    break
            if not is_dup:
                all_verts.append(v)
    return all_verts
def analyze_face_geometry(face, num_segments=CURVE_SEGMENTS):
    alibre_area = None
    try:
        alibre_area = face.GetArea()
        if alibre_area == 0:
            alibre_area = None
    except Exception:
        pass
    edges = []
    try:
        edges = face.GetEdges() or []
    except Exception:
        pass
    face_vertices = []
    try:
        face_vertices = face.GetVertices() or []
    except Exception:
        pass
    is_rect = False
    try:
        is_rect = face.IsRectangle()
    except Exception:
        pass
    edge_infos = [EdgeInfo(e, i) for i, e in enumerate(edges)]
    circular_edges = [e for e in edge_infos if e.is_circular]
    curved_edges = [e for e in edge_infos if e.is_curved and not e.is_circular]
    linear_edges = [e for e in edge_infos if not e.is_curved and not e.is_circular and e.vertex_count >= 2]
    diameters = [e.diameter for e in circular_edges if e.diameter]
    unique_diameters = list(set([round(d, 6) for d in diameters]))
    if is_rect and len(linear_edges) == 4:
        lengths = sorted([e.length for e in linear_edges if e.length])
        if len(lengths) == 4:
            width = (lengths[0] + lengths[1]) / 2.0
            height = (lengths[2] + lengths[3]) / 2.0
            return ("rectangle", RectangleProperties(width, height), alibre_area)
    if len(circular_edges) == 1 and is_complete_circle(circular_edges[0]) and len(linear_edges) == 0:
        radius = circular_edges[0].diameter / 2.0
        return ("circle", CircleProperties(radius), alibre_area)
    if len(circular_edges) == 2 and all(is_complete_circle(e) for e in circular_edges):
        r1 = circular_edges[0].diameter / 2.0
        r2 = circular_edges[1].diameter / 2.0
        outer_r, inner_r = max(r1, r2), min(r1, r2)
        if outer_r - inner_r < 0.001:
            return ("circle", CircleProperties(outer_r), alibre_area)
        return ("annulus", AnnulusProperties(outer_r, inner_r), alibre_area)
    if len(unique_diameters) == 1 and len(circular_edges) >= 1 and len(linear_edges) == 0:
        radius = unique_diameters[0] / 2.0
        total_arc = sum(e.length for e in circular_edges if e.length)
        expected_circumference = math.pi * unique_diameters[0]
        if abs(total_arc - expected_circumference) < expected_circumference * 0.02:
            return ("circle", CircleProperties(radius), alibre_area)
    if len(unique_diameters) == 2:
        r1, r2 = max(unique_diameters) / 2.0, min(unique_diameters) / 2.0
        if r1 - r2 < 0.001:
            return ("circle", CircleProperties(r1), alibre_area)
        return ("annulus", AnnulusProperties(r1, r2), alibre_area)
    if len(circular_edges) == 1 and len(linear_edges) == 1:
        arc_edge = circular_edges[0]
        if arc_edge.diameter and arc_edge.length:
            expected_semicircle = math.pi * arc_edge.diameter / 2.0
            if abs(arc_edge.length - expected_semicircle) < expected_semicircle * ARC_LENGTH_TOLERANCE:
                radius = arc_edge.diameter / 2.0
                return ("semicircle", SemicircleProperties(radius), alibre_area)
    if len(circular_edges) == 1 and len(linear_edges) == 2:
        arc_edge = circular_edges[0]
        if arc_edge.diameter and arc_edge.length:
            expected_quarter = math.pi * arc_edge.diameter / 4.0
            if abs(arc_edge.length - expected_quarter) < expected_quarter * ARC_LENGTH_TOLERANCE:
                radius = arc_edge.diameter / 2.0
                return ("quarter_circle", QuarterCircleProperties(radius), alibre_area)
    if len(circular_edges) == 2 and len(linear_edges) == 2:
        d1 = circular_edges[0].diameter
        d2 = circular_edges[1].diameter
        if d1 and d2 and abs(d1 - d2) < 0.01:
            l1 = circular_edges[0].length
            l2 = circular_edges[1].length
            expected_semi = math.pi * d1 / 2.0
            if l1 and l2:
                if abs(l1 - expected_semi) < expected_semi * ARC_LENGTH_TOLERANCE and abs(l2 - expected_semi) < expected_semi * ARC_LENGTH_TOLERANCE:
                    lengths = sorted([e.length for e in linear_edges if e.length])
                    if len(lengths) == 2 and abs(lengths[0] - lengths[1]) < 0.01:
                        slot_length = lengths[0] + d1
                        slot_width = d1
                        return ("slot", SlotProperties(slot_length, slot_width), alibre_area)
    if len(linear_edges) == 8 and len(circular_edges) == 0:
        loops = find_edge_loops(linear_edges)
        if len(loops) == 2:
            dims1 = get_rectangle_dimensions(loops[0])
            dims2 = get_rectangle_dimensions(loops[1])
            if dims1 and dims2:
                area1 = dims1[0] * dims1[1]
                area2 = dims2[0] * dims2[1]
                if area1 > area2:
                    outer_w, outer_h = dims1
                    inner_w, inner_h = dims2
                else:
                    outer_w, outer_h = dims2
                    inner_w, inner_h = dims1
                return ("rectangular_tubing", RectangularTubingProperties(outer_w, outer_h, inner_w, inner_h), alibre_area)
    if len(linear_edges) >= 6 and len(circular_edges) == 0:
        loops = find_edge_loops(linear_edges)
        if len(loops) == 2 and len(loops[0]) >= 3 and len(loops[1]) >= 3:
            outer_verts = collect_unique_vertices(loops[0])
            inner_verts = collect_unique_vertices(loops[1])
            if len(outer_verts) >= 3 and len(inner_verts) >= 3:
                outer_2d = project_to_2d(outer_verts)
                inner_2d = project_to_2d(inner_verts)
                outer_poly = PolygonProperties(outer_2d, "outer")
                inner_poly = PolygonProperties(inner_2d, "inner")
                if outer_poly.area > inner_poly.area:
                    shape_name = "Polygonal Tube (%d/%d sides)" % (len(loops[0]), len(loops[1]))
                    return ("polygonal_tube", PolygonWithPolygonalHoleProperties(outer_2d, inner_2d, shape_name), alibre_area)
                else:
                    shape_name = "Polygonal Tube (%d/%d sides)" % (len(loops[1]), len(loops[0]))
                    return ("polygonal_tube", PolygonWithPolygonalHoleProperties(inner_2d, outer_2d, shape_name), alibre_area)
        elif len(loops) >= 3:
            loop_data = []
            for loop in loops:
                if len(loop) >= 3:
                    verts = collect_unique_vertices(loop)
                    if len(verts) >= 3:
                        try:
                            verts_2d = project_to_2d(verts)
                            poly = PolygonProperties(verts_2d, "temp")
                            loop_data.append((poly.area, verts_2d))
                        except Exception:
                            pass
            if len(loop_data) >= 2:
                loop_data.sort(key=lambda x: x[0], reverse=True)
                outer_2d = loop_data[0][1]
                inner_boundaries = [ld[1] for ld in loop_data[1:]]
                shape_name = "Ship Section (%d openings)" % len(inner_boundaries)
                return ("ship_section", ShipSectionProperties(outer_2d, inner_boundaries, shape_name=shape_name), alibre_area)
    all_edges = linear_edges + curved_edges + circular_edges
    if len(all_edges) >= 6:
        loops = find_edge_loops(all_edges)
        if len(loops) >= 3:
            loop_data = []
            for loop in loops:
                boundary_pts = []
                for edge_info in loop:
                    if edge_info.is_curved or edge_info.is_circular:
                        pts = discretize_curved_edge(edge_info, max(8, num_segments // 8))
                        for pt in pts:
                            is_dup = any(vec_distance(pt, ex) < VERTEX_TOLERANCE for ex in boundary_pts)
                            if not is_dup:
                                boundary_pts.append(pt)
                    else:
                        for v in edge_info.vertices:
                            is_dup = any(vec_distance(v, ex) < VERTEX_TOLERANCE for ex in boundary_pts)
                            if not is_dup:
                                boundary_pts.append(v)
                if len(boundary_pts) >= 3:
                    try:
                        verts_2d = project_to_2d(boundary_pts)
                        poly = PolygonProperties(verts_2d, "temp")
                        loop_data.append((poly.area, verts_2d))
                    except Exception:
                        pass
            if len(loop_data) >= 2:
                loop_data.sort(key=lambda x: x[0], reverse=True)
                outer_2d = loop_data[0][1]
                inner_boundaries = [ld[1] for ld in loop_data[1:]]
                shape_name = "Ship Section (%d openings, curved)" % len(inner_boundaries)
                return ("ship_section", ShipSectionProperties(outer_2d, inner_boundaries, shape_name=shape_name), alibre_area)
    if len(linear_edges) >= 3 and len(circular_edges) == 1 and is_complete_circle(circular_edges[0]):
        hole_edge = circular_edges[0]
        if hole_edge.diameter:
            hole_radius = hole_edge.diameter / 2.0
            outer_verts = collect_unique_vertices(linear_edges)
            if len(outer_verts) >= 3:
                outer_2d = project_to_2d(outer_verts)
                shape_name = "Polygon (%d sides) with Hole" % len(linear_edges)
                return ("polygon_with_hole", PolygonWithCircularHoleProperties(outer_2d, hole_radius, shape_name=shape_name), alibre_area)
    if len(linear_edges) == 3 and len(circular_edges) == 0:
        all_verts = collect_unique_vertices(linear_edges)
        if len(all_verts) == 3:
            points_2d = project_to_2d(all_verts)
            return ("triangle", PolygonProperties(points_2d, "Triangle"), alibre_area)
    if len(linear_edges) >= 3 and len(circular_edges) == 0 and len(curved_edges) == 0:
        all_verts = collect_unique_vertices(linear_edges)
        if len(all_verts) >= 3:
            points_2d = project_to_2d(all_verts)
            return ("polygon", PolygonProperties(points_2d, "Polygon (%d sides)" % len(linear_edges)), alibre_area)
    if len(curved_edges) > 0 or (len(circular_edges) > 0 and not all(is_complete_circle(e) for e in circular_edges)):
        boundary_points = discretize_face_boundary(edge_infos, num_segments)
        if len(boundary_points) >= 3:
            points_2d = project_to_2d(boundary_points)
            return ("freeform", PolygonProperties(points_2d, "Freeform (%d pts)" % len(points_2d)), alibre_area)
    boundary_points = []
    for v in face_vertices:
        try:
            pt = [v.X, v.Y, v.Z]
            is_dup = any(vec_distance(pt, ex) < VERTEX_TOLERANCE for ex in boundary_points)
            if not is_dup:
                boundary_points.append(pt)
        except Exception:
            pass
    if len(boundary_points) >= 3:
        points_2d = project_to_2d(boundary_points)
        return ("polygon", PolygonProperties(points_2d, "Polygon (from vertices)"), alibre_area)
    if alibre_area and alibre_area > 0:
        radius = math.sqrt(alibre_area / math.pi)
        props = CircleProperties(radius)
        props.shape_type = "Equivalent Circle (from area)"
        return ("equivalent_circle", props, alibre_area)
    return ("unknown", None, alibre_area)
def discretize_face_boundary(edge_infos, num_segments):
    boundary_points = []
    ordered_edges = order_edges_by_connectivity(edge_infos)
    for edge_info in ordered_edges:
        if is_complete_circle(edge_info) and edge_info.diameter:
            radius = edge_info.diameter / 2.0
            center = estimate_circle_center(edge_info)
            for i in range(num_segments):
                angle = 2.0 * math.pi * i / num_segments
                pt = [center[0] + radius * math.cos(angle), center[1] + radius * math.sin(angle), center[2]]
                boundary_points.append(pt)
        elif edge_info.is_curved or edge_info.is_circular:
            pts = discretize_curved_edge(edge_info, max(8, num_segments // 4))
            for pt in pts:
                is_dup = any(vec_distance(pt, ex) < VERTEX_TOLERANCE for ex in boundary_points)
                if not is_dup:
                    boundary_points.append(pt)
        else:
            for v in edge_info.vertices:
                is_dup = any(vec_distance(v, ex) < VERTEX_TOLERANCE for ex in boundary_points)
                if not is_dup:
                    boundary_points.append(v)
    return boundary_points
def order_edges_by_connectivity(edge_infos):
    if not edge_infos:
        return []
    ordered = []
    remaining = list(edge_infos)
    ordered.append(remaining.pop(0))
    while remaining:
        last_edge = ordered[-1]
        found_idx = -1
        for i, edge in enumerate(remaining):
            if edges_share_vertex(last_edge, edge):
                found_idx = i
                break
        if found_idx >= 0:
            ordered.append(remaining.pop(found_idx))
        else:
            ordered.append(remaining.pop(0))
    return ordered
def edges_share_vertex(edge1, edge2):
    for v1 in edge1.vertices:
        for v2 in edge2.vertices:
            if vec_distance(v1, v2) < VERTEX_TOLERANCE:
                return True
    return False
def estimate_circle_center(edge_info):
    if edge_info.vertices:
        n = len(edge_info.vertices)
        cx = sum(v[0] for v in edge_info.vertices) / n
        cy = sum(v[1] for v in edge_info.vertices) / n
        cz = sum(v[2] for v in edge_info.vertices) / n
        return [cx, cy, cz]
    return [0.0, 0.0, 0.0]
def discretize_curved_edge(edge_info, num_points):
    points = []
    num_points = max(8, num_points)
    if edge_info.is_circular and edge_info.diameter and len(edge_info.vertices) >= 2:
        v1, v2 = edge_info.vertices[0], edge_info.vertices[1]
        chord = vec_distance(v1, v2)
        radius = edge_info.diameter / 2.0
        if chord < 2.0 * radius - 1e-9:
            mid_chord = vec_midpoint(v1, v2)
            chord_vec = vec_subtract(v2, v1)
            chord_len = vec_length(chord_vec)
            if chord_len > 1e-9:
                perp_2d = [-chord_vec[1], chord_vec[0], 0.0]
                perp_len = vec_length(perp_2d)
                if perp_len > 1e-9:
                    perp_2d = vec_scale(perp_2d, 1.0 / perp_len)
                    sagitta = radius - math.sqrt(max(0, radius**2 - (chord/2.0)**2))
                    arc_mid = vec_add(mid_chord, vec_scale(perp_2d, sagitta))
                    center = vec_add(mid_chord, vec_scale(perp_2d, -(radius - sagitta)))
                    v1_rel = vec_subtract(v1, center)
                    v2_rel = vec_subtract(v2, center)
                    angle1 = math.atan2(v1_rel[1], v1_rel[0])
                    angle2 = math.atan2(v2_rel[1], v2_rel[0])
                    if angle2 < angle1:
                        angle2 += 2.0 * math.pi
                    arc_angle = angle2 - angle1
                    if arc_angle > math.pi:
                        angle1, angle2 = angle2, angle1 + 2.0 * math.pi
                        arc_angle = angle2 - angle1
                    for i in range(num_points + 1):
                        t = i / float(num_points)
                        angle = angle1 + t * arc_angle
                        pt = [center[0] + radius * math.cos(angle),
                              center[1] + radius * math.sin(angle),
                              center[2]]
                        points.append(pt)
                    return points
    if len(edge_info.vertices) >= 2:
        v1, v2 = edge_info.vertices[0], edge_info.vertices[1]
        if edge_info.is_curved and edge_info.length and edge_info.chord_length > 1e-9:
            bulge_ratio = (edge_info.length / edge_info.chord_length - 1.0) * 0.5
            mid = vec_midpoint(v1, v2)
            chord_vec = vec_subtract(v2, v1)
            perp = [-chord_vec[1], chord_vec[0], chord_vec[2]]
            perp_len = vec_length(perp)
            if perp_len > 1e-9:
                perp = vec_scale(perp, bulge_ratio * edge_info.chord_length / perp_len)
                control = vec_add(mid, perp)
                for i in range(num_points + 1):
                    t = i / float(num_points)
                    p01 = vec_lerp(v1, control, t)
                    p12 = vec_lerp(control, v2, t)
                    pt = vec_lerp(p01, p12, t)
                    points.append(pt)
                return points
        for i in range(num_points + 1):
            t = i / float(num_points)
            points.append(vec_lerp(v1, v2, t))
    return points
def project_to_2d(points_3d):
    n = len(points_3d)
    if n < 3:
        raise ValueError("Need at least 3 points for 2D projection")
    cx = sum(p[0] for p in points_3d) / n
    cy = sum(p[1] for p in points_3d) / n
    cz = sum(p[2] for p in points_3d) / n
    origin = [cx, cy, cz]
    p0 = points_3d[0]
    v1 = None
    for i in range(1, n):
        vec = vec_subtract(points_3d[i], p0)
        if vec_length(vec) > VERTEX_TOLERANCE:
            v1 = vec
            break
    if v1 is None:
        raise ValueError("All points are coincident")
    v2 = None
    for i in range(2, n):
        vec = vec_subtract(points_3d[i], p0)
        cross = vec_cross(v1, vec)
        if vec_length(cross) > VERTEX_TOLERANCE:
            v2 = vec
            break
    if v2 is None:
        if abs(v1[2]) < 0.9:
            v2 = vec_cross(v1, [0.0, 0.0, 1.0])
        else:
            v2 = vec_cross(v1, [1.0, 0.0, 0.0])
    normal = vec_normalize(vec_cross(v1, v2))
    u_axis = vec_normalize(v1)
    v_axis = vec_cross(normal, u_axis)
    points_2d = []
    for p3d in points_3d:
        rel = vec_subtract(p3d, origin)
        u = vec_dot(rel, u_axis)
        v = vec_dot(rel, v_axis)
        points_2d.append([u, v])
    cx2 = sum(p[0] for p in points_2d) / n
    cy2 = sum(p[1] for p in points_2d) / n
    def angle_key(p):
        return math.atan2(p[1] - cy2, p[0] - cx2)
    points_2d.sort(key=angle_key)
    return points_2d
def fmt(value, decimals=6):
    if value is None:
        return "N/A"
    if isinstance(value, float) and (math.isinf(value) or math.isnan(value)):
        return "N/A"
    if abs(value) < 1e-10:
        return "0.0"
    elif abs(value) >= 1e6 or (abs(value) < 0.001 and abs(value) > 1e-10):
        return "%.4e" % value
    else:
        return ("%%.%df" % decimals) % value
UNIT_CONVERSIONS = {
    'mm': 1.0,
    'cm': 0.1,
    'm': 0.001,
    'in': 1.0 / 25.4,
    'ft': 1.0 / 304.8
}
def convert_value(value, power, unit):
    if value is None or (isinstance(value, float) and (math.isinf(value) or math.isnan(value))):
        return value
    factor = UNIT_CONVERSIONS[unit] ** power
    return value * factor
def generate_face_section(face_name, props, alibre_area, unit):
    u1 = unit
    u2 = unit + "^2"
    u3 = unit + "^3"
    u4 = unit + "^4"
    p = props
    lines = []
    lines.append("  %s (%s)" % (face_name, props.shape_type))
    lines.append("  " + "-" * 40)
    if hasattr(props, 'radius') and not hasattr(props, 'outer_radius'):
        lines.append("    Radius:      %s %s" % (fmt(convert_value(props.radius, 1, unit)), u1))
    if hasattr(props, 'outer_radius'):
        lines.append("    Outer R:     %s %s" % (fmt(convert_value(props.outer_radius, 1, unit)), u1))
        lines.append("    Inner R:     %s %s" % (fmt(convert_value(props.inner_radius, 1, unit)), u1))
    if hasattr(props, 'width') and hasattr(props, 'height'):
        lines.append("    Width:       %s %s" % (fmt(convert_value(props.width, 1, unit)), u1))
        lines.append("    Height:      %s %s" % (fmt(convert_value(props.height, 1, unit)), u1))
    if hasattr(props, 'outer_width') and hasattr(props, 'outer_height'):
        lines.append("    Outer W:     %s %s" % (fmt(convert_value(props.outer_width, 1, unit)), u1))
        lines.append("    Outer H:     %s %s" % (fmt(convert_value(props.outer_height, 1, unit)), u1))
    if hasattr(props, 'inner_width') and hasattr(props, 'inner_height'):
        lines.append("    Inner W:     %s %s" % (fmt(convert_value(props.inner_width, 1, unit)), u1))
        lines.append("    Inner H:     %s %s" % (fmt(convert_value(props.inner_height, 1, unit)), u1))
    if hasattr(props, 'hole_radius'):
        lines.append("    Hole R:      %s %s" % (fmt(convert_value(props.hole_radius, 1, unit)), u1))
    if hasattr(props, 'semi_major') and hasattr(props, 'semi_minor'):
        lines.append("    Semi-Major:  %s %s" % (fmt(convert_value(props.semi_major, 1, unit)), u1))
        lines.append("    Semi-Minor:  %s %s" % (fmt(convert_value(props.semi_minor, 1, unit)), u1))
    if hasattr(props, 'outer_major') and hasattr(props, 'inner_major'):
        lines.append("    Outer Major: %s %s" % (fmt(convert_value(props.outer_major, 1, unit)), u1))
        lines.append("    Outer Minor: %s %s" % (fmt(convert_value(props.outer_minor, 1, unit)), u1))
        lines.append("    Inner Major: %s %s" % (fmt(convert_value(props.inner_major, 1, unit)), u1))
        lines.append("    Inner Minor: %s %s" % (fmt(convert_value(props.inner_minor, 1, unit)), u1))
    if hasattr(props, 'length') and hasattr(props, 'width') and not hasattr(props, 'height'):
        lines.append("    Length:      %s %s" % (fmt(convert_value(props.length, 1, unit)), u1))
        lines.append("    Width:       %s %s" % (fmt(convert_value(props.width, 1, unit)), u1))
    if hasattr(props, 'flange_width') and hasattr(props, 'web_thickness'):
        lines.append("    Height:      %s %s" % (fmt(convert_value(props.height, 1, unit)), u1))
        lines.append("    Flange W:    %s %s" % (fmt(convert_value(props.flange_width, 1, unit)), u1))
        lines.append("    Web t:       %s %s" % (fmt(convert_value(props.web_thickness, 1, unit)), u1))
        lines.append("    Flange t:    %s %s" % (fmt(convert_value(props.flange_thickness, 1, unit)), u1))
    if hasattr(props, 'leg_height') and hasattr(props, 'leg_width'):
        lines.append("    Leg H:       %s %s" % (fmt(convert_value(props.leg_height, 1, unit)), u1))
        lines.append("    Leg W:       %s %s" % (fmt(convert_value(props.leg_width, 1, unit)), u1))
        lines.append("    Thickness:   %s %s" % (fmt(convert_value(props.thickness, 1, unit)), u1))
    if hasattr(props, 'gap'):
        lines.append("    Gap:         %s %s" % (fmt(convert_value(props.gap, 1, unit)), u1))
    if hasattr(props, 'top_width') and hasattr(props, 'bottom_width'):
        lines.append("    Top W:       %s %s" % (fmt(convert_value(props.top_width, 1, unit)), u1))
        lines.append("    Bottom W:    %s %s" % (fmt(convert_value(props.bottom_width, 1, unit)), u1))
    if hasattr(props, 'lip_height'):
        lines.append("    Lip H:       %s %s" % (fmt(convert_value(props.lip_height, 1, unit)), u1))
    if hasattr(props, 'top_flange_width'):
        lines.append("    Top Flg W:   %s %s" % (fmt(convert_value(props.top_flange_width, 1, unit)), u1))
        lines.append("    Top Flg t:   %s %s" % (fmt(convert_value(props.top_flange_thickness, 1, unit)), u1))
        lines.append("    Bot Flg W:   %s %s" % (fmt(convert_value(props.bottom_flange_width, 1, unit)), u1))
        lines.append("    Bot Flg t:   %s %s" % (fmt(convert_value(props.bottom_flange_thickness, 1, unit)), u1))
    if hasattr(props, 'top_thickness') and hasattr(props, 'bottom_thickness'):
        lines.append("    Top t:       %s %s" % (fmt(convert_value(props.top_thickness, 1, unit)), u1))
        lines.append("    Bottom t:    %s %s" % (fmt(convert_value(props.bottom_thickness, 1, unit)), u1))
        lines.append("    Left t:      %s %s" % (fmt(convert_value(props.left_thickness, 1, unit)), u1))
        lines.append("    Right t:     %s %s" % (fmt(convert_value(props.right_thickness, 1, unit)), u1))
    if hasattr(props, 'hole_diameter'):
        lines.append("    Hole D:      %s %s" % (fmt(convert_value(props.hole_diameter, 1, unit)), u1))
        lines.append("    Spacing X:   %s %s" % (fmt(convert_value(props.hole_spacing_x, 1, unit)), u1))
        lines.append("    Spacing Y:   %s %s" % (fmt(convert_value(props.hole_spacing_y, 1, unit)), u1))
        lines.append("    Num Holes:   %d" % props.num_holes)
    if hasattr(props, 'num_openings'):
        lines.append("    Openings:    %d" % props.num_openings)
    if hasattr(props, 'neutral_axis_height'):
        lines.append("    NA Height:   %s %s" % (fmt(convert_value(props.neutral_axis_height, 1, unit)), u1))
    if hasattr(props, 'c_deck') and hasattr(props, 'c_keel'):
        lines.append("    c (deck):    %s %s" % (fmt(convert_value(props.c_deck, 1, unit)), u1))
        lines.append("    c (keel):    %s %s" % (fmt(convert_value(props.c_keel, 1, unit)), u1))
    lines.append("    Area:        %s %s" % (fmt(convert_value(p.area, 2, unit)), u2))
    lines.append("    Ix-x:        %s %s" % (fmt(convert_value(p.Ix_centroid, 4, unit)), u4))
    lines.append("    Iy-y:        %s %s" % (fmt(convert_value(p.Iy_centroid, 4, unit)), u4))
    lines.append("    J:           %s %s" % (fmt(convert_value(p.J_centroid, 4, unit)), u4))
    lines.append("    rx-x:        %s %s" % (fmt(convert_value(p.rx, 1, unit)), u1))
    lines.append("    ry-y:        %s %s" % (fmt(convert_value(p.ry, 1, unit)), u1))
    lines.append("    Sx-x (min):  %s %s" % (fmt(convert_value(p.Sx_min, 3, unit)), u3))
    lines.append("    Sy-y (min):  %s %s" % (fmt(convert_value(p.Sy_min, 3, unit)), u3))
    if hasattr(p, 'Sx_deck') and hasattr(p, 'Sx_keel'):
        lines.append("    Sx (deck):   %s %s" % (fmt(convert_value(p.Sx_deck, 3, unit)), u3))
        lines.append("    Sx (keel):   %s %s" % (fmt(convert_value(p.Sx_keel, 3, unit)), u3))
    lines.append("")
    return "\n".join(lines)
def generate_batch_report(face_results):
    lines = []
    lines.append("")
    lines.append("=" * 70)
    lines.append("            BATCH AREA MOMENTS REPORT")
    lines.append("            %s v%s" % (ScriptName, ScriptVersion))
    lines.append("=" * 70)
    lines.append("")
    lines.append("Faces Analyzed: %d" % len(face_results))
    lines.append("")
    for unit in ['mm', 'cm', 'm', 'in', 'ft']:
        lines.append("=" * 70)
        lines.append("RESULTS IN %s" % unit.upper())
        lines.append("=" * 70)
        lines.append("")
        for face_name, props, alibre_area in face_results:
            if props is not None:
                lines.append(generate_face_section(face_name, props, alibre_area, unit))
    lines.append("=" * 70)
    lines.append("END OF REPORT")
    lines.append("=" * 70)
    return "\n".join(lines)
def get_all_faces(part):
    faces = []
    face_index = 0
    try:
        all_faces = part.GetFaces()
        if all_faces:
            for face in all_faces:
                face_index += 1
                try:
                    name = None
                    if hasattr(face, 'Name'):
                        name = face.Name
                    if not name:
                        name = "Face<%d>" % face_index
                    faces.append((name, face))
                except Exception:
                    faces.append(("Face<%d>" % face_index, face))
    except Exception:
        pass
    def sort_key(x):
        name = x[0]
        if name.startswith("Face<") and name.endswith(">"):
            try:
                return (0, int(name[5:-1]))
            except ValueError:
                return (1, name)
        return (1, name)
    faces.sort(key=sort_key)
    return faces
def run():
    try:
        part = CurrentPart()
        if part is None:
            print "ERROR: Please open a part before running this script."
            return
    except Exception:
        print "ERROR: Please open a part before running this script."
        return
    all_faces = get_all_faces(part)
    if not all_faces:
        print "ERROR: No faces found in the part. Make sure the part contains solid geometry."
        return
    print "=" * 70
    print "BATCH AREA MOMENTS - Scanning %d faces..." % len(all_faces)
    print "=" * 70
    face_results = []
    for face_name, face in all_faces:
        print "Analyzing: %s" % face_name
        try:
            shape_type, props, alibre_area = analyze_face_geometry(face, CURVE_SEGMENTS)
            if props is not None:
                if shape_type not in ("polygon_with_hole", "circle", "annulus", "rectangle", "rectangular_tubing",
                                      "ellipse", "elliptical_tube", "quarter_circle", "slot", "polygonal_tube",
                                      "i_section", "t_section", "channel", "angle", "z_section",
                                      "hat", "sigma", "omega", "double_angle", "double_channel",
                                      "built_up_i", "built_up_box", "perforated_plate", "composite", "ship_section"):
                    if alibre_area and alibre_area > 0 and props.area > 0:
                        area_ratio = alibre_area / props.area
                        if area_ratio > 1.10 or area_ratio < 0.90:
                            print "  -> Skipping (non-planar, area mismatch: %.1f%%)" % ((area_ratio - 1.0) * 100)
                            continue
                face_results.append((face_name, props, alibre_area))
                print "  -> %s (Area: %.2f mm^2)" % (props.shape_type, props.area)
            else:
                print "  -> Could not determine geometry"
        except Exception as ex:
            print "  -> Error: %s" % str(ex)
    print ""
    if not face_results:
        print "ERROR: No faces could be analyzed. Make sure faces are planar."
        return
    report = generate_batch_report(face_results)
    print report
    sys.stdout.flush()
    print ""
    print "=" * 70
    print "Batch analysis complete! %d faces analyzed." % len(face_results)
    print "=" * 70
run()

@stephensmitchell

Copy link
Copy Markdown
Author

This is a debug script, so the length isn't as important, but it shouldn't grow more than this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment