Developing a generalized cone.

[edwilliams16]

There's ongoing interest in creating a developable connection between two curves - eg a sheet metal fitting between two arbitrary pipe sections.

A developable surface is one that can be flattened without any stretching or compression. It is one that be created by suitably bending a flat sheet. The two common examples are the curved surfaces of cones and cylinders. The most general developable surface is a union of portions of generalized cylinders and generalized cones.
A generalized cylinder is the extrusion of a space curve in a fixed direction. For a standard cylinder, the curve is a circle (or circular arc).
Similarly, a generalized cone is a space curve lofted to a point with a ruled surface.

All developable surfaces are ruled surfaces - but not all ruled surfaces are developable. The simplest counter-example is the ruled surface loft between two non-parallel line segments.

In viewtopic.php?p=178312#p178312 @Chris_G showed how to flatten a generalized cylinder.
His idea was to re-parameterize the surface with u distance around the cross-section and v distance along axial direction. Mapping the edges onto a plane unrolls the surface.

For a generalized cone, I re-parameterized the surface in polar coordinates. u being angle and v radius. Then the flattened surface is plotted on the XY-plane by converting Polar to Cartesian. I think my code can be simplified - but for now, here it is:

Code: Select all

#start with edge and vertex
#create loft to point from edge to vertex
#intersect it with sphere from vertex
#reparameterize with length, then angle as u coordinate
#create an expanded loft to point,  overlapping edge
#reparametrize v to length
#find u, v coords of edge on the surface and draw flattened triangle in 2d.
from math import sin, cos
debug = False
doc = App.ActiveDocument()
#input is an edge an a vertex
def flattenGeneralizedCone(edge, vertex):
    #edge = doc.getObject('Sketch').Shape.Edge1
    #vertex = doc.getObject('Vertex').Shape.Vertex1
    if edge.Curve.TypeId == 'Part::GeomBSplineCurve':
        skspline = edge.Curve
    else:
        skspline = edge.Curve.toBSpline(*edge.ParameterRange)

    vskspline = skspline.copy()
    for i, pole in enumerate(vskspline.getPoles()):
        vskspline.setPole(i+1, vertex.Point) #note 1-origin for pole list

    #create regular ruled surface from skspline to vertex
    sksurf = Part.BSplineSurface()
    sksurf.buildFromNSections([vskspline, skspline])
    if debug: Part.show(sksurf.toShape(), 'sksurf')

    #we are going to make a new surface that includes the sketch
    radii = [(v - vertex.Point).Length for v in skspline.discretize(20)]
    minrad = min(radii)
    maxrad = max(radii)

    radius = minrad/10
    center = vertex.Point
    sph = Part.makeSphere(radius, center)
    sphsurf = sph.Face1.Surface
    sphcurve = sphsurf.intersect(sksurf)[0]# this curve is constant radius from vertex on sksurf

    #create a version of sphcurve parameterized by length along the curve
    npts = 20
    dis = sphcurve.length()/npts
    pts = sphcurve.discretize(Distance = dis)
    params = [i * dis for i in range(len(pts))]
    sphcurvei = Part.BSplineCurve()
    sphcurvei.interpolate(Points = pts, Parameters = params)
    #now we create a surface coincident with the loft with polar coordinates u =r v=theta
    vsphcurvei = sphcurvei.copy()
    for i, pole in enumerate(sphcurvei.getPoles()):
        vsphcurvei.setPole(i+1, vertex.Point)

    osphcurvei = sphcurvei.copy()
    for i, pole in enumerate(sphcurvei.getPoles()):
        v = vertex.Point + (pole - vertex.Point) * (maxrad/minrad) *20
        osphcurvei.setPole(i+1, v)

    newsurf = Part.BSplineSurface()
    newsurf.buildFromNSections([vsphcurvei, osphcurvei])
    newsurfshp = newsurf.toShape()
    newsurf.scaleKnotsToBounds(0, newsurfshp.Edge4.Length/newsurfshp.Edge1.Length, 0, newsurfshp.Edge1.Length) #polar coords u is theta v is r
    newsurfshp = newsurf.toShape()
    if debug: 
        newsurfdoc = Part.show(newsurfshp, 'newsurf')
        newsurfdoc.ViewObject.Visibility = False

    #map skspline to the u,v plane
    skpts = skspline.discretize(Distance = 0.1)
    skpcoords = [newsurf.parameter(p) for p in skpts]
    skccoords = [App.Vector(uv[1]*cos(uv[0]), uv[1] * sin(uv[0]), 0.) for uv in skpcoords]
    flatbs = Part.BSplineCurve()
    flatbs.interpolate(skccoords)
    #close up in the uv plane
    l1 = Part.LineSegment(App.Vector(), flatbs.StartPoint).toShape()
    l2 = Part.LineSegment(flatbs.EndPoint, App.Vector()).toShape()
    wire = Part.Wire([l1, flatbs.toShape(), l2])
    Part.show(Part.Face(wire), "Flattened Face")

    #fit a 2d spline to the uv figure
    bsonuv = Part.Geom2d.BSplineCurve2d()
    bsonuv.interpolate([App.Base.Vector2d(uv[0], uv[1]) for uv in skpcoords])
    edgeonnewsurf = bsonuv.toShape(newsurf, bsonuv.FirstParameter, bsonuv.LastParameter)
    if debug: Part.show(edgeonnewsurf, "edgeonnewsurf")

    l3 = Part.Geom2d.Line2dSegment(App.Base.Vector2d(skpcoords[0][0],0), App.Base.Vector2d(*skpcoords[0]))
    l4 = Part.Geom2d.Line2dSegment(App.Base.Vector2d(*skpcoords[-1]), App.Base.Vector2d(0,0))
    l3s = l3.toShape(newsurf, l3.FirstParameter, l3.LastParameter)
    l4s = l4.toShape(newsurf, l4.FirstParameter, l4.LastParameter)
    wireonface = Part.Wire(Part.__sortEdges__([l3s, l4s,edgeonnewsurf]))
    if debug: Part.show(wireonface, 'WireOnFace')

    #map back onto newsurf
    #newsurf is coincident with the original loft, but has polar coordinate parameters.
    ltp = Part.Face(newsurf, wireonface)
    ltp.validate()
    Part.show(ltp, "LoftToPoint")


edgecount = 0
vertexcount = 0
for sel in Gui.Selection.getSelectionEx('',0):
    for path in sel.SubElementNames if sel.SubElementNames else ['']:
        obj = sel.Object.getSubObject(path)
        if obj.ShapeType == 'Edge':
            edge = obj
            edgecount += 1
        elif obj.ShapeType == 'Vertex':
            vertex = obj
            vertexcount += 1

if edgecount != 1 or vertexcount != 1 or len(edge.Vertexes) != 2:
    App.Console.PrintMessage("select one open edge and one vertex\n")
else:
    flattenGeneralizedCone(edge, vertex)

Select an open edge and a vertex and a developable surface will be created, with polar coordinate parameterization. The flattened surface is created on the XY-Plane. I need to special case closed edges.

With this one has a necessary ingredient to triangulate the loft between two arbitrary space -curves, making a developable surface.
However, the flat pieces would need to be glued together, and some automated way of choosing the triangulation would be needed. I have some ideas...

To use this in order to add capability to Curves|Flatten Face we need a simple way of identifying surfaces as generalized cones or cylinders. Sampling the Gaussian Curvature is one way, but I think given a ruled surface, a 1d sampling along an edge may suffice. The intersection of two lines a principal curvature direction could find the vertex.

Screenshot 2023-03-14 at 4.55.56 PM.png (20.37 KiB) Viewed 1010 times

[jfc4120]

I also have a general cone program (symmetrical only). It just lays out a regular cone, thought I put it here if anyone needs it:

Code: Select all

import FreeCAD, FreeCADGui
# -*- coding: utf-8 -*-
import FreeCAD, FreeCADGui
import Draft
import math
# from math import cos, sin, radians
from PySide import QtGui

convert = 25.4
reduct =  .039370078
TOTAL = 0
getDouble = QtGui.QInputDialog.getDouble
Vector, Placement = App.Vector, App.Placement
doc = App.ActiveDocument

start_point = Vector(0, 0, 0)
group = App.ActiveDocument.addObject("App::DocumentObjectGroupPython", "Group")


def line_length(x1=0, y1=0, z1=0, length=10, angle=0):
    x2 = x1 + length * math.cos(math.radians(angle))
    y2 = y1 + length * math.sin(math.radians(angle))
    z2 = z1
    #Draft.makeWire([Vector(x1, y1, z1), Vector(x2, y2, z2)])
    group.addObjects([Draft.makeWire([Vector(x1, y1, z1), Vector(x2, y2, z2)])])

def coord(x1=0, y1=0, z1=0, length=10, angle=0):
    x2 = x1 + length * math.cos(math.radians(angle))
    x2 = x2 * reduct
    y2 = y1 + length * math.sin(math.radians(angle))
    y2 = y2 * reduct
    z2 = z1
    #z2 = z2 * convert
    coo = [x2, y2, z2]
    return coo

D, dflag = getDouble(None, 'example', 'Large diameter', decimals=3)
C, cflag = getDouble(None, 'example', 'Small diameter', decimals=3)
B, bflag = getDouble(None, 'example', 'Height', decimals=3)
#D = 20  #base
#C = 10 #top
#B = 15 #height

R = math.sqrt((D * B / (D - C)) ** 2 + (D / 2) ** 2)
print(R)

RP = R * C / D

THETA = 180 * D / R
print(THETA)

zaxis = App.Vector(0, 0, 1)
p1 = App.Vector(0, 0, 0)
place1 = App.Placement(p1, App.Rotation(zaxis, 0))
circle1 = Draft.make_circle(R * convert, placement=place1, face=None, startangle=0, endangle=THETA, support=None)

p2 = App.Vector(0, 0, 0)
place2 = App.Placement(p1, App.Rotation(zaxis, 0))
circle2 = Draft.make_circle(RP * convert, placement=place1, face=None, startangle=0, endangle=THETA, support=None)

x1 = y1 = z1 = 0  # Edit coordinate origin
length = D
length = length * 25.4
angle = 180
line_length(x1, y1, z1, length, angle)
clist1 = coord(x1, y1, z1, length, angle)

x1 = (-1 * (D / 2)) + (C / 2)
x1 = x1 * convert

y1 = B * convert
z1 = 0  # Edit coordinate origin
length = C
length = length * 25.4
angle = 180
line_length(x1, y1, z1, length, angle)
clist2 = coord(x1, y1, z1, length, angle)

group.addObjects([Draft.makeWire([Vector(0, 0, 0), Vector(x1, y1, z1)])])

group.addObjects([Draft.makeWire([Vector(-D * convert, 0, 0), Vector(clist2[0] * convert, clist2[1] * convert, clist2[2] * convert)])])
ZP = R - RP

hline = Draft.make_line(Vector(RP * convert, 0, 0), Vector(R * convert, 0, 0))

X = ((RP * math.cos(math.radians(THETA))))
Y = ((RP * math.sin(math.radians(THETA))))
X1 = ((R * math.cos(math.radians(THETA))))
Y1 = ((R * math.sin(math.radians(THETA))))

hline = Draft.make_line(Vector(X * convert , Y * convert, 0), Vector(X1 * convert, Y1 * convert, 0))

Probably needs some cleanup, so tweak as needed.

[edwilliams16]

I assume this is for a regular circle cross section cone? My code will work with an arbitrary BSpline.
@jfc4120

[jfc4120]

@edwilliams16 yes it just lays out a cone that is symmetrical.

But also @nm2107 on his post viewtopic.php?t=76594 has also made some breakthroughs on his program.
My program depends on setting points but it will layout (unfold) any round to round with any miter. But is is still in work. I am still converting from basiccad to python.

Also does your program require the curves workbemch?

I generally like straight analytic geometry and trigonometry solutions, that way any layout program can easily be switched over to another language.

In other words I generally do fittings using basic distance between points and angle between lines.

I admit that I have studied your square to round program, and it is still way over my head. But as I keep at it, it is starting to fall into place.

Edit:

Another program I have will layout the first piece of a round elbow, or select 2 piece you can get a pipe miter, looks like:

In inches, but someone can convert to mm in the code.

Also in your cone program, is a small slit required to unfold like you need to unwrap a mesh?

I only included this here since it's also round type work, and if downloading yours they may want to download cylindrical also.

[Chris_G]

In viewtopic.php?p=178312#p178312 @Chris_G showed how to flatten a generalized cylinder.
His idea was to re-parameterize the surface with u distance around the cross-section and v distance along axial direction. Mapping the edges onto a plane unrolls the surface.

I did experiment with flattening surfaces of extrusion a couple of weeks ago.
I used the same workflow as yours.
Unfortunately, I lack free time nowadays to push this further.
But here is my latest code.
It still fails on seam edges of closed faces.
Select a face of type SurfaceOfExtrusion and run.

Code: Select all

def arclength_approx(curve, num=100):
    "Returns an arc-length approximation BSpline of a curve or edge"
    sh = curve
    if hasattr(curve, "toShape"):
        sh = curve.toShape()
    dist = sh.Length / num
    pts = sh.discretize(Distance=dist)
    params = [i * dist for i in range(len(pts) - 1)]
    params.append(sh.Length)
    bs = Part.BSplineCurve()
    bs.approximate(Points=pts, Parameters=params)
    return bs

# get the face of an extrusion
s = FreeCADGui.Selection.getSelectionEx()
face = s[0].SubObjects[0]

u0, u1, v0, v1 = face.ParameterRange
surf = face.Surface
if hasattr(surf, "BasisSurface"):
    surf = surf.BasisSurface

basc = surf.BasisCurve
eps = 0  #1e-2
based = basc.toShape(u0 - eps, u1 + eps)
# Part.show(based, "BasisCurve")

# create a normal plane face
pl = Part.Plane(basc.value(basc.FirstParameter), surf.Direction)
plts = Part.RectangularTrimmedSurface(pl, -1e3, 1e3, -1e3, 1e3)
plsh = plts.toShape()
# Part.show(plsh)

# project basis curve on normal plane
proj = plsh.project([based])  # , surf.Direction)
# Part.show(proj, "BasisCurve Projection")

if len(proj.Edges) > 1:
    w = Part.Wire(proj.Edges)
    pe = w.approximate(1e-10, 1e-7, 1000, 7)
    print("Projection : multiple edges")
else:
    pe = proj.Edge1.Curve

# create a arc-lenght approx of projected curve
alc = arclength_approx(pe)
# Part.show(alc.toShape(), "BasisCurve Approx")

# extrude arc-lenght parametrized curve
sof = Part.SurfaceOfExtrusion(alc, surf.Direction)
u0, u1, v0, v1 = sof.bounds()
nts = Part.RectangularTrimmedSurface(sof, u0, u1, -1e3, 1e3)
nface = nts.toShape()
# Part.show(nface)

# Transfer original face edges on XY plane
xy = Part.Plane()
fl = []
for w in face.Wires:
    for e in w.Edges:
        pr = nface.project([e])
        if len(pr.Edges) > 0:
            cos, fp, lp = nface.curveOnSurface(pr.Edge1)
            print(cos, fp, lp)
            fl.append(cos.toShape(xy, fp, lp))

# try to build a face or a wire
print(len(fl))
flsh = Part.Wire(fl)
if flsh.isValid():
    flfa = Part.Face(flsh)
    flfa.validate()
    if flfa.isValid():
        Part.show(flfa, "Flat Face")
    else:
        Part.show(flsh, "Flat Wire")
else:
    comp = Part.Compound(fl)
    Part.show(comp, "Flat Compound")

[edwilliams16]

I did experiment with flattening surfaces of extrusion a couple of weeks ago.

Any ideas how best to identify a generalized cone/loft to a point? TheirTypeId is 'Part::GeomBSplineSurface', which isn't helpful.

[edwilliams16]

I did experiment with flattening surfaces of extrusion a couple of weeks ago.

Any ideas how best to identify a generalized cone/loft to a point? TheirTypeId is 'Part::GeomBSplineSurface', which isn't helpful.


I've made some progress in creating developable lofts - so far restricted to open convex curves on parallel planes.

Screenshot 2023-03-15 at 12.23.34 PM.png (23.93 KiB) Viewed 773 times

[onekk]

...

Maybe not relevant bu I remeber an old post of wmayer, were he show a way to know if a point is comprised in a surface, using UV coordinates, using as example a Sphere.

There were a method to check if a uv coordinate is in the "surface limits".

Another thing maybe related to the above I have seen that if a point is very near the boundaries, probably due to the approximation, sometimes it fails to project wires on a surface, noted this near the seam line of a Cylinder surface, I don't know if it relevant here, it will be interesting to investigate this behaviour, sadly, I have no code around that could be used as MWE.

Sorry if above is not relevant.

Kind Regards

Carlo D.

[edwilliams16]

Screenshot 2023-03-16 at 12.38.37 PM.png (13.01 KiB) Viewed 636 times

But here is my latest code.
It still fails on seam edges of closed faces.
Select a face of type SurfaceOfExtrusion and run.

I worked on your code a little. I think I've fixed the co-planar type missing edge problems by extrapolating the approx curve. I've yet to look at extruded loops and seams.

Code: Select all

def arclength_approx(curve, num=100, extend = [True, True]):
    '''Returns an arc-length approximation BSpline of a curve or edge
       if extend Extrapolate ends by an extra point'''
    sh = curve
    if hasattr(curve, "toShape"):
        sh = curve.toShape()
    dist = sh.Length / num
    if debug:
        print([v.Point for v in sh.Vertexes])

    deltau = abs((sh.LastParameter - sh.FirstParameter)/ num ) 
    # not quite right - not equal distance
    if extend[0]:
        shd = sh.Curve.getD3(sh.FirstParameter)
        leftpt = [shd[0] - deltau * (shd[1] - deltau/2*(shd[2] + deltau/3 * shd[3]))] #cubic Taylor series
        if debug:
            Part.show(Part.Vertex(leftpt[0]), "leftpt")
            #print(leftpt[0])
    else:
        leftpt = []
    if extend[1]:
        shd = sh.Curve.getD3(sh.LastParameter)
        rightpt = [shd[0] + deltau * (shd[1] + deltau/2*(shd[2] + deltau/3 * shd[3]))]
        if debug:
            Part.show(Part.Vertex(rightpt[0]), 'rightpt')
            #print(rightpt[0])
    else:
        rightpt = []

    pts = leftpt 
    pts += sh.discretize(Distance=dist)
    pts += rightpt
    params = [i * dist for i in range(len(pts))]
    bs = Part.BSplineCurve()
    bs.approximate(Points=pts, Parameters=params)
    return bs


debug = True
# get the face of an extrusion
s = FreeCADGui.Selection.getSelectionEx()
face = s[0].SubObjects[0]

u0, u1, v0, v1 = face.ParameterRange
surf = face.Surface
if hasattr(surf, "BasisSurface"):
    surf = surf.BasisSurface

basc = surf.BasisCurve
eps = 0
based = basc.toShape(u0 - eps, u1 + eps)
if debug: Part.show(based, "BasisCurve")

# create a normal plane face
pl = Part.Plane(basc.value(basc.FirstParameter), surf.Direction)
plts = Part.RectangularTrimmedSurface(pl, -1e3, 1e3, -1e3, 1e3)
plsh = plts.toShape()
#if debug: Part.show(plsh)

# project basis curve on normal plane
proj = plsh.project([based])  # , surf.Direction)
if debug: Part.show(proj, "BasisCurve Projection")

if len(proj.Edges) > 1:
    w = Part.Wire(proj.Edges)
    pe = w.approximate(1e-10, 1e-7, 1000, 7)
    print("Projection : multiple edges")
else:
    pe = proj.Edge1.Curve

# create a arc-lenght approx of projected curve
alc = arclength_approx(pe)
if debug: Part.show(alc.toShape(), "BasisCurve Approx")

# extrude arc-lenght parametrized curve
sof = Part.SurfaceOfExtrusion(alc, surf.Direction)
u0, u1, v0, v1 = sof.bounds()
nts = Part.RectangularTrimmedSurface(sof, u0, u1, -1e2, 1e2)
nface = nts.toShape()
if debug: Part.show(nface)

# Transfer original face edges on XY plane
xy = Part.Plane()
fl = []
for w in face.Wires:
    for e in w.Edges:
        pr = nface.project([e])
        if len(pr.Edges) > 0:
            cos, fp, lp = nface.curveOnSurface(pr.Edge1)
            print(cos, fp, lp)
            fl.append(cos.toShape(xy, fp, lp))

# try to build a face or a wire
#print(len(fl))
flsortedges = Part.sortEdges(fl)
flwires = [Part.Wire(s) for s in flsortedges]
flface = Part.Face(flwires)
flface.validate()
if flface.isValid():
    Part.show(flface, "FlatFace")
else:
    comp = Part.Compound(flwires)
    Part.show(comp, FlatWires)

Screenshot 2023-03-16 at 12.38.37 PM.png (13.01 KiB) Viewed 636 times