This was interesting to me, so I poked around, adding some diagnostics to the code.
Code: Select all
"""Sample code.
This code was written as an sample code
Name: filename.py
Author: Carlo Dormeletti
Copyright: 2022
Licence: CC BY-NC-ND 4.0 IT
"""
import os # noqa
import sys # noqa
import math # noqa
from math import cos, pi, sin, sqrt # noqa
import FreeCAD
import FreeCADGui # noqa
from FreeCAD import Placement, Rotation, Vector # noqa
import Part # noqa
DOC = FreeCAD.activeDocument()
DOC_NAME = "test_offset"
def clear_doc():
"""Clear active document deleting all the objects."""
for obj in DOC.Objects:
DOC.removeObject(obj.Name)
def setview():
"""Rearrange View."""
FreeCAD.Gui.SendMsgToActiveView("ViewFit")
FreeCAD.Gui.activeDocument().activeView().viewAxometric()
if DOC is None:
FreeCAD.newDocument(DOC_NAME)
FreeCAD.setActiveDocument(DOC_NAME)
DOC = FreeCAD.activeDocument()
else:
clear_doc()
# EPS= tolerance to use to cut the parts
EPS = 0.10
EPS_C = EPS * -0.5
ROT0 = Rotation(0, 0, 0)
# CODE START HERE
def tangent2C(c1, c2, rad1, rad2):
"""Calculate four tangent point to join two circles.
Courtesy of edwilliams16
Parameters:
c1 Vector Center of circle1
c2 Vector Center of circle2
rad1 float Radius of circle1
rad2 float Radius of circle2
Return:
[tg1, tg2, tg3, tg4] tangent points.
If one circle strictly encloses the other or are they are concentric
return []
"""
lg = (c1 - c2).Length
if abs(rad1 - rad2) > lg or lg == 0:
return []
else:
c = (rad1 - rad2) / lg
s = sqrt(1 - c * c)
n = (c2 - c1) / lg # unit vector along line of centers
m = Vector(0, 0, 1).cross(n) # orthogonal vector
vp = c * n + s * m # unit vectors from center to tangent points
vm = c * n - s * m
return [c1 + rad1 * vp, c1 + rad1 * vm, c2 + rad2 * vp, c2 + rad2 * vm]
def dest_pt(pt, angle, length):
"""Calculate destination point.
Vector math courtesy of edwilliams16
Parameters:
pt = starting point
angle = rad
length = units
"""
# return Vector(pt.x + math.cos(angle) * length, pt.y + math.sin(angle) * length, pt.z)
lv = Vector(length, 0, 0)
p_rot = Rotation(math.degrees(angle), 0, 0)
dept = Placement(pt, p_rot).multVec(lv)
return dept
def mk_shell():
"""Make an shell for test."""
c1cen = Vector(0, 0, 0)
c1rad = 100
c2cen = Vector(0, 100, 0)
c2rad = 50
tgp1, tgp2, tgp3, tgp4 = tangent2C(c1cen, c2cen, c1rad, c2rad)
cir1 = Part.makeCircle(c1rad, c1cen)
cir2 = Part.makeCircle(c2rad, c2cen)
edg1 = Part.LineSegment(tgp3, tgp1).toShape()
edg2 = Part.LineSegment(tgp2, tgp4).toShape()
# obtain cir1 section between tangents
par1 = cir1.Curve.parameter(tgp1)
par2 = cir1.Curve.parameter(tgp2)
tgv1 = cir1.Curve.value(par1)
tgv2 = cir1.Curve.value(par2)
c1_3pos = tgv1 - (tgv2 - tgv1) * 0.5
par3 = cir1.Curve.parameter(c1_3pos)
tgv3 = cir1.Curve.value(par3)
edg3 = Part.ArcOfCircle(tgv2, tgv3, tgv1).toShape()
# obtain cir2 section between tangents
par4 = cir2.Curve.parameter(tgp3)
par5 = cir2.Curve.parameter(tgp4)
tgv4 = cir2.Curve.value(par4)
tgv5 = cir2.Curve.value(par5)
c2_3pos = tgv4 + (tgv5 - tgv4) * 0.5
par6 = cir2.Curve.parameter(c2_3pos)
tgv6 = cir2.Curve.value(par6)
edg4 = Part.ArcOfCircle(tgv4, tgv6, tgv5).toShape()
# edge are ordered in counterclockwise
obj_f = Part.Wire((edg2, edg3, edg1, edg4))
# Part.show(obj_f, "wire")
return obj_f
def find_inters():
"""Find intersection between lines and a shell."""
n_sp = 21
ang_adv = 360 / n_sp
ln_lg = 200
sh_thi = 80
s_lines = []
for idx in range(0, n_sp):
end_pt = dest_pt(
Vector(0, 0, 0),
pi / 2.0 + math.radians(ang_adv * idx),
ln_lg)
line = Part.LineSegment(
Vector(0, 0, sh_thi * 0.5),
Vector(end_pt.x, end_pt.y, sh_thi * 0.5))
s_lines.append(line)
Part.show(line.toShape(), "line{}".format(idx))
wire = mk_shell()
shell = wire.extrude(Vector(0, 0, sh_thi))
Part.show(shell, "shell")
diag_len = shell.BoundBox.DiagonalLength
edges = []
for line in s_lines:
edges.append(line.toShape(0.0, diag_len))
comp = Part.Compound(edges)
d, ipts, info = comp.distToShape(shell)
#print(f'd[0]= {d[0]} ipts[0]= {ipts[0]} info = {info[0]}')
points_on_lines, points_on_shell = list(zip(*ipts))
Part.show(Part.Compound([Part.Vertex(p) for p in points_on_lines]))
Part.show(Part.Compound([Part.Vertex(p) for p in points_on_shell]))
return (d, ipts, info, comp)
d, ipts, info, comp = find_inters()
setview()
#look at lineN intersection with the shell
N = 5
lN = App.ActiveDocument.getObject('line'+str(N)).Shape
shl= App.ActiveDocument.getObject('shell').Shape
shp=App.ActiveDocument.getObject('Shape').Shape
infoN = info[N]
print(f'Intersection from Shape list at {shp.Vertexes[N].Point}')
print(f'Intersection on line at {lN.valueAt(infoN[2])}')
face = shl.Faces[infoN[4]]
faceuv = infoN[5] #u,v coordinate of point on the shell face
print(f'Intersection on shell Face at {face.valueAt(*faceuv)}')
edge = comp.Edges[infoN[1]]
edgeu = infoN[2] # u coordinate along edge
print(f'Intersection on edge at {edge.valueAt(edgeu)}')
# in addition to valueAt() we have normalAt() froviding a normal vector
# and derivative1At() providing tangent vectors
The info return locates the intersections (Actually nearest points) on the two shapes - in this case lines and a shell, giving the coordinate(s) on the edge or face.
The sample code shows the numerous ways of extracting the location of the point. By the same route, you can find the normal and tangent vectors at the point, which you can use to construct further geometry.
There were two puzzles that came up:
(1) d was returned a single scalar. Shouldn't is be a list of distances for the items in the compound?
(2) In the case of the face, derivative2At() provides (I think) d2
x/du^2 and d2
x/dv^2. If this is so, the cross-derivative d2
x/dudv is missing, and would be required to characterize the surface curvature.
I've not followed up on these yet.