.DXF files

[samhughes-cfd]

Hi all,

I am using this excellent package for all sort of section property calculations, mainly focusing on geometric properties of hydrofoils.

The limiting condition on the accuracy of this 2D FEM code for me is the preparation of the .dxf which is to be meshed. I'm wanting to ask on open question how people prepare their .dxf files for this program?

Currently I'm extracting projected sketches of cross sections from Abaqus; exporting as .stp. A very robust procedure. Then I am manually tracing a spline curve in SolidEdge to then export this curve as a .dxf in the save in flat option from the metal sheet work environment. The point density is up to the clicks of my cursor and eye-judgement from trialling different point density patterns. This part is awfully slow and unreliable... some of my .dxfs are processed perfectly, some are processed quite erroneously with a larger delta area between the Abaqus computed value of area mm^2 vs. sectionproperties computed value of area mm^2.

Does anyone know a robust and automated way to get the .dxfs from a .stp? Abaqus also exports sketch as .igs and .sat.

Finally, I appreciate this is a slightly adjacent issue to the sectionproperties package. But prepration of .dxf for me is the limiting factor on my accuracy right now!

Example of erroneous .dxf:

Example of non erroneous .dxf:

Best regards all!

[github-actions]

Thanks for opening your first issue in sectionproperties 🙌 Pull requests are always welcome 😉

[samhughes-cfd]

I should share the files I am using to plot the raw geometry, then process into polylines. Apologies the code blocks arent working too well.

section_calc_n/utils/raw_geometry.py, this is what is visualised in the multicoloured plot

`import logging
import ezdxf
import numpy as np
from pathlib import Path
from typing import List, Tuple

class RawGeometry:
def init(
self,
filepath: str | Path,
label: str,
*,
n_line: int = 50,
n_poly: int = 10,
logs_dir: Path = None
) -> None:
self.filepath = Path(filepath)
self.label = label
self.n_line = max(n_line, 2)
self.n_poly = max(n_poly, 2)
self._raw_pts: List[List[Tuple[float, float]]] = []

    # ───── Set up logger ─────
    self.logger = logging.getLogger(f"RawGeometry.{label}")
    self.logger.setLevel(logging.INFO)

    if logs_dir:
        logs_dir.mkdir(parents=True, exist_ok=True)
        log_path = logs_dir / f"raw_geometry_{label}.log"
        file_handler = logging.FileHandler(log_path, mode="w", encoding="utf-8")
        formatter = logging.Formatter("%(asctime)s - %(levelname)s - %(message)s")
        file_handler.setFormatter(formatter)
        if not self.logger.handlers:
            self.logger.addHandler(file_handler)

def _densify(self, p0, p1, n) -> List[Tuple[float, float]]:
    xs = np.linspace(p0[0], p1[0], n)
    ys = np.linspace(p0[1], p1[1], n)
    return list(zip(xs, ys))

def extract(self) -> List[List[Tuple[float, float]]]:
    if self._raw_pts:
        return self._raw_pts

    self.logger.info("Reading DXF file: %s", self.filepath)
    try:
        doc = ezdxf.readfile(self.filepath)
    except Exception as e:
        self.logger.error("Failed to read DXF file %s: %s", self.filepath, e)
        raise

    entity_count = 0
    for ent in doc.modelspace():
        typ = ent.dxftype()
        if typ not in {"LINE", "LWPOLYLINE", "POLYLINE", "SPLINE", "ARC"}:
            continue

        try:
            if typ == "LINE":
                p0 = (ent.dxf.start.x, ent.dxf.start.y)
                p1 = (ent.dxf.end.x, ent.dxf.end.y)
                pts = self._densify(p0, p1, self.n_line)

            elif typ in {"LWPOLYLINE", "POLYLINE"}:
                verts = [(x, y) for x, y, *_ in ent.get_points()]
                pts = []
                for a, b in zip(verts, verts[1:]):
                    seg_pts = self._densify(a, b, self.n_poly)[:-1]
                    pts.extend(seg_pts)
                pts.append(verts[-1])

            elif typ == "SPLINE":
                pts = [(x, y) for x, y, *_ in ent.approximate(200)]

            elif typ == "ARC":
                th = np.linspace(np.radians(ent.dxf.start_angle), np.radians(ent.dxf.end_angle), 100)
                cx, cy, r = ent.dxf.center.x, ent.dxf.center.y, ent.dxf.radius
                pts = [(cx + r * np.cos(t), cy + r * np.sin(t)) for t in th]

            self._raw_pts.append(pts)
            entity_count += 1

        except Exception as exc:
            self.logger.debug("Skipped malformed entity %s: %s", typ, exc)

    if not self._raw_pts:
        raise ValueError(f"No drawable entities in {self.filepath}")

    self.logger.info("Raw geometry extraction complete. %d drawable entities parsed.", entity_count)
    return self._raw_pts`

section_calc_n/utils/processed_geometry.py, this is what is viewed in the plot with the black dots. Note errors arise from this step depending on the .DXF characteristics and this propagates into the meshed domain and thus my answers

`from pathlib import Path
import logging
import numpy as np
from shapely.geometry import Polygon
from shapely.affinity import scale
from sectionproperties.pre.geometry import Geometry, CompoundGeometry
from geometry_transforms.twist_offset import TwistOffset
from geometry_transforms.centroid_offset import CentroidOffset

def _cosine_resample_exterior(poly: Polygon, n_total: int) -> Polygon:
if n_total < 3:
raise ValueError("Need at least 3 exterior points")

xy = np.asarray(poly.exterior.coords[:-1])
seg_len = np.linalg.norm(np.diff(xy, axis=0), axis=1)
cum_len = np.concatenate([[0.0], np.cumsum(seg_len)])
length = cum_len[-1]
t = 0.5 * (1 - np.cos(np.linspace(0.0, np.pi, n_total)))
s_targets = t * length

new_ext = []
j = 0
for s in s_targets:
    while s > cum_len[j + 1]:
        j += 1
    frac = (s - cum_len[j]) / seg_len[j]
    pt = (1 - frac) * xy[j] + frac * xy[j + 1]
    new_ext.append(tuple(pt))

return Polygon(new_ext, [ring.coords[:] for ring in poly.interiors])

class ProcessedGeometry:
def init(
self,
filepath: Path,
label: str,
logs_dir: Path,
*,
spline_delta: float = 0.05,
degrees_per_segment: float = 0.5,
exterior_nodes: int = 400,
):
self.filepath = Path(filepath)
self.label = label
self.logs_dir = logs_dir
self.spline_delta = spline_delta
self.degrees_per_segment = degrees_per_segment
self.exterior_nodes = exterior_nodes
self.geometry: Geometry | CompoundGeometry | None = None

def extract_and_transform(self, twist_deg: float, cx: float, cy: float) -> Geometry | CompoundGeometry:
    logging.info("[%s] Importing DXF: '%s'", self.label, self.filepath)

    geom_raw = Geometry.from_dxf(
        dxf_filepath=self.filepath,
        spline_delta=self.spline_delta,
        degrees_per_segment=self.degrees_per_segment
    )

    if isinstance(geom_raw, Geometry): 
        scaled_geom = scale(geom_raw.geom, xfact=0.001, yfact=0.001, origin=(0, 0)) # .dxf in mm scale to m
        geom_raw = Geometry(geom=scaled_geom)
    elif isinstance(geom_raw, CompoundGeometry): 
        scaled_geoms = []
        for g in geom_raw.geoms:
            scaled = scale(g.geom, xfact=0.001, yfact=0.001, origin=(0, 0)) # .dxf in mm scale to m
            scaled_geoms.append(Geometry(geom=scaled))
        geom_raw = CompoundGeometry(scaled_geoms)
    else:
        raise TypeError("Unexpected geometry type when scaling.")

    raw_polys = []
    if isinstance(geom_raw, Geometry):
        polys = geom_raw.geom.geoms if geom_raw.geom.geom_type == "MultiPolygon" else [geom_raw.geom]
        raw_polys.extend(polys)
    else:
        for g in geom_raw.geoms:
            polys = g.geom.geoms if g.geom.geom_type == "MultiPolygon" else [g.geom]
            raw_polys.extend(polys)

    dense_polys = []
    for i, p in enumerate(raw_polys):
        try:
            dense_polys.append(_cosine_resample_exterior(p, self.exterior_nodes))
        except Exception as e:
            logging.warning("[%s] Resampling polygon %d failed: %s", self.label, i, e)

    geom_list = [Geometry(geom=p) for p in dense_polys]
    geom = geom_list[0] if len(geom_list) == 1 else CompoundGeometry(geom_list)

    # ───── Rotate ───────────────────────────────────────────── COMMENT IF NO TWIST
    geom_rotated = TwistOffset(
        geometry=geom,
        desired_twist_deg=twist_deg,
        label=self.label,
        logs_dir=self.logs_dir,
        twist_axis_ratio= 0.333  # centre of twist at 0.33*chord  from leading edge
    ).apply()


    # ───── Translate ──────────────────────────────────────────
    geom_translated = CentroidOffset(
        geometry=geom_rotated,    # with no twist comment replace "geometry=geom_rotated" with "geometry=geom"
        cx_target=cx,
        cy_target=cy,
        label=self.label,
        logs_dir=self.logs_dir
    ).apply()

    self.geometry = geom_translated
    return self.geometry`