solver, singularities, and macros

[dlajoie]

Hello,

I am trying to optimize lengths on a figure. A sample file is provided, along with a short screencast to illustrate the issue.

http://lutolf.vserver.nimag.net/files/t ... o-show.mp4
http://lutolf.vserver.nimag.net/files/tri-macro.fcstd

The problem is quite straight forward: consider an equilateral triangle of which the COG coincides with the COG of an arbitrary triangle of fixed position. The equilateral triangle can rotate, in the plane, around its COG. The corners of the arbitrary triangles are linked to the center of the sides of the equilateral triangle, one by one, using lines whose lenghts need to be optimized : that is, the total length of these 3 lines must be at a minimum.


There are many problems associated with this : first of all, setting a datum on one of these lines "constrains" the sketch (according to the solver, because in fact we know in this case the sketch is not fully constrained: there are usually two different rotations of the equilateral triangle for a given line length (unless the datum is the max or min allowable value in which case there is only one valid position, or unless the datum is beyond the min/max allowable values in which case there are none valid positions)

This leads straight to the root of the problem : trying to set a datum for any other of the white segments returns the error "cannot set negative value for datum XX" (XX being - probably? - the id of the previous datum set). This happens every time, even when the default (measured) value is not changed.

This is a pure case of "not-so-redundant constraints": while there are two valid sketch positions with a single datum set, there will be only two valid values for the second constraint/datum, either of which _will_ effectively constrain the sketch.



I am trying to write a macro to optimize those lengths, and this issue is really a problem. Since I cannot set the datum for the two remaining lines, I also cannot _read_ their values (trying to set the datum does show me the actual value (length) of each line, but the script doesn't _get_ this value and I _really_ don't want having to set it by hand).



Another representation of the underlying problem is seen at the start of the video : trying to make many-turns rotation of the equilateral triangle shows the links (white lines) do not quite behave like one should expect they would - although the sketch is supposedly fully constrained (except for the rotation of the equi-triangle), the solver (?) runs into problems because of the symmetry of the figure.


It looks like the solver obviously does not go "beyond singularity" when checking for constrained sketches. As a consequence, I suppose FreeCAD does not use quaternions for plane (flat) figures (probably not 3D figures either). This might solve that.


Since I'm not very optimistic about the real issue being fixed anytime soon, suggestions for getting around the symptoms faced with the macro are particularly welcome.

[microelly2]

You can yuse SketchObjectPython instead of a simple SketchObject.
then you "move your points" near the expected solution behind the singularity and run the solve method of the sketch.
some times it helps to add/delete helper constraints on the fly.

[TheMarkster]

The macro sets an angle constraint on line 27 (index 26) in a loop from 0 degrees to 360 degrees in 1 degree increments. You can change the precision if you want to narrow down the minimum. You should change the range of start / stop angles, too, since by then you'll know the generally where the optimal is at (somewhere around 85 degrees, give or take, with a total length of 106.938 for all 3 lines).

The gif was recorded at 2 frames per second to keep it under the 1 mb limit for attaching on here. The chart was generated with Open Office spreadsheet.

tri-macro.gif (732.63 KiB) Viewed 610 times

tri-macro.png (57.51 KiB) Viewed 610 times

Code: Select all

import FreeCAD
import math
from PySide import QtGui
from time import sleep

def processEvents():
     sleep(0.001)
     QtGui.QApplication.processEvents()
def getDistance3d(x1, y1, z1, x2, y2, z2):
    return math.sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)

def dist(p1,p2):
    return getDistance3d(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z)

def lineLength(line):
    return dist(line.StartPoint,line.EndPoint)

sketch = App.ActiveDocument.Sketch005
geom = sketch.Geometry
con = sketch.Constraints

topLine = geom[30]
rightLine = geom[24]
leftLine = geom[23]

topLineLength = lineLength(topLine)
rightLineLength = lineLength(rightLine)
leftLineLength = lineLength(leftLine)
totalLength = topLineLength+rightLineLength+leftLineLength

#sketch.setDatum(76,App.Units.Quantity('12.500000 deg'))

minAngle = 0
maxAngle = 360

#minAngle = 30
#maxAngle = 35

maxTotal= 0.0
maxTotalAngle = minAngle
minTotal = 100000
minTotalAngle = minAngle
precision = 1.0 #change to 10.0, 100.0, etc. for better precision (but will take longer)

csvString=[]

for angle in range(int(minAngle*precision),int(maxAngle*precision),1):
    sketch = App.ActiveDocument.Sketch005
    geom = sketch.Geometry
    con = sketch.Constraints

    topLine = geom[30]
    rightLine = geom[24]
    leftLine = geom[23]
    try:
        fAngle = float(angle/precision)#fAngle = float(angle/10.0) for more precision
        sketch.addConstraint(Sketcher.Constraint('Angle',26,fAngle*math.pi/180.0))
        sketch.setDatum(76,App.Units.Quantity(str(fAngle)+' deg'))
    except:
        FreeCAD.Console.PrintMessage("Exception with angle = "+str(fAngle)+' deg\n')
        sketch.delConstraint(76)
    App.ActiveDocument.recompute()
    processEvents()
    topLineLength = lineLength(topLine)
    rightLineLength = lineLength(rightLine)
    leftLineLength = lineLength(leftLine)
    totalLength = topLineLength+rightLineLength+leftLineLength
    if totalLength > maxTotal:
        maxTotal = totalLength
        maxTotalAngle = fAngle
        FreeCAD.Console.PrintMessage("MaxTotal: "+str(maxTotal)+" was found with angle = "+str(float(maxTotalAngle))+"\n")
    if totalLength < minTotal:
        minTotal = totalLength
        minTotalAngle = fAngle
        FreeCAD.Console.PrintMessage("MinTotal: "+str(minTotal)+" was found with angle = "+str(float(minTotalAngle))+"\n")
    FreeCAD.Console.PrintMessage("angle = "+str(fAngle)+"\n"+"top line = "+str(topLineLength)+"\n right line = "+str(rightLineLength)+"\n left line = "+str(leftLineLength)+"\n total = "+str(totalLength)+"\n")
    csvString.append(str(fAngle)+','+str(totalLength)+','+str(topLineLength)+','+str(rightLineLength)+','+str(leftLineLength)+'\n')

    try:
        sketch.delConstraint(76)
    except:
        FreeCAD.Console.PrintMessage("Exception when deleting constraint 76\n")
#FreeCAD.Console.PrintMessage("top line = "+str(topLineLength)+"\n right line = "+str(rightLineLength)+"\n left line = "+str(leftLineLength)+"\n total = "+str(totalLength)+"\n")

    FreeCAD.Console.PrintMessage("MaxTotal: "+str(maxTotal)+" was found with angle = "+str(float(maxTotalAngle))+"\n")
    FreeCAD.Console.PrintMessage("MinTotal: "+str(minTotal)+" was found with angle = "+str(float(minTotalAngle))+"\n")

    FreeCAD.Console.PrintMessage("CSV:\n\n")
    FreeCAD.Console.PrintMessage("Angle,Total Length,Top Line, Right Line, Left Line\n")
    for csv in csvString:
        FreeCAD.Console.PrintMessage(csv)

[dlajoie]

Hey!

TheMarkster: thanks for the macro, this seems to be it. I will try in the next few days (maybe even hours?)

The graph definitely shows what is going on. I was not expecting the individual line lengths to be in phase, though.


microelly2: I didn't really take the time to test your solution. it wasn't totally clear for me and I had to put the project on hold anyway, so I can't really say how valid it was.



either way, I will read TheMarkster's code to see how it works.

I totally thank you for your help and hope you had fun thinking about or solving the problem. Whenever the time is right (I hope as soon as possible) I will tell you more about why this is. It may be that you like the project


Happy hacking



ps: I am quite surprised how different the behaviour is with altering the figure with a) the mouse b) the script. If anyone can enlighten me on why it is so smooth with the script and so messed up with the mouse, I am curious the read why,

[TheMarkster]

Has been a while since I worked on that. I've learned a few things since then. Could no doubt have been done more robustly. Seems it doesn't work with 0.18 because the new version finds a redundant constraint that 0.17 didn't catch. Consequently, you have to delete the redundant constraint or else the solver doesn't update the sketch, and no animation. But since I used direct indexing (index 76) to refer to the new angle constraint, the macro fails. Quick fix is to delete the redundant constraint and use 75 in place of 76.

Code: Select all

import FreeCAD
import math
from PySide import QtGui
from time import sleep

def processEvents():
     sleep(0.001)
     #QtGui.qApp.processEvents() #so FreeCAD stays responsive to user input
     QtGui.QApplication.processEvents()
def getDistance3d(x1, y1, z1, x2, y2, z2):
    return math.sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)

def dist(p1,p2):
    return getDistance3d(p1.x,p1.y,p1.z,p2.x,p2.y,p2.z)

def lineLength(line):
    return dist(line.StartPoint,line.EndPoint)

sketch = App.ActiveDocument.Sketch005
geom = sketch.Geometry
con = sketch.Constraints

topLine = geom[30]
rightLine = geom[24]
leftLine = geom[23]

topLineLength = lineLength(topLine)
rightLineLength = lineLength(rightLine)
leftLineLength = lineLength(leftLine)
totalLength = topLineLength+rightLineLength+leftLineLength

#sketch.setDatum(76,App.Units.Quantity('12.500000 deg'))

minAngle = 0
maxAngle = 360

#minAngle = 30
#maxAngle = 35

maxTotal= 0.0
maxTotalAngle = minAngle
minTotal = 100000
minTotalAngle = minAngle
precision = 1.0 #change to 10.0, 100.0, etc. for better precision (but will take longer)

csvString=[]
constraint=75
for angle in range(int(minAngle*precision),int(maxAngle*precision),1):
    sketch = App.ActiveDocument.Sketch005
    geom = sketch.Geometry
    con = sketch.Constraints

    topLine = geom[30]
    rightLine = geom[24]
    leftLine = geom[23]
    try:
        fAngle = float(angle/precision)#fAngle = float(angle/10.0) for more precision
        sketch.addConstraint(Sketcher.Constraint('Angle',26,fAngle*math.pi/180.0))
        sketch.setDatum(constraint,App.Units.Quantity(str(fAngle)+' deg'))
    except:
        FreeCAD.Console.PrintMessage("Exception with angle = "+str(fAngle)+' deg\n')
        sketch.delConstraint(constraint)
    App.ActiveDocument.recompute()
    processEvents()
    topLineLength = lineLength(topLine)
    rightLineLength = lineLength(rightLine)
    leftLineLength = lineLength(leftLine)
    totalLength = topLineLength+rightLineLength+leftLineLength
    if totalLength > maxTotal:
        maxTotal = totalLength
        maxTotalAngle = fAngle
        FreeCAD.Console.PrintMessage("MaxTotal: "+str(maxTotal)+" was found with angle = "+str(float(maxTotalAngle))+"\n")
    if totalLength < minTotal:
        minTotal = totalLength
        minTotalAngle = fAngle
        FreeCAD.Console.PrintMessage("MinTotal: "+str(minTotal)+" was found with angle = "+str(float(minTotalAngle))+"\n")
    FreeCAD.Console.PrintMessage("angle = "+str(fAngle)+"\n"+"top line = "+str(topLineLength)+"\n right line = "+str(rightLineLength)+"\n left line = "+str(leftLineLength)+"\n total = "+str(totalLength)+"\n")
    csvString.append(str(fAngle)+','+str(totalLength)+','+str(topLineLength)+','+str(rightLineLength)+','+str(leftLineLength)+'\n')

    try:
        sketch.delConstraint(constraint)
    except:
        FreeCAD.Console.PrintMessage("Exception when deleting constraint 76\n")
#FreeCAD.Console.PrintMessage("top line = "+str(topLineLength)+"\n right line = "+str(rightLineLength)+"\n left line = "+str(leftLineLength)+"\n total = "+str(totalLength)+"\n")

    FreeCAD.Console.PrintMessage("MaxTotal: "+str(maxTotal)+" was found with angle = "+str(float(maxTotalAngle))+"\n")
    FreeCAD.Console.PrintMessage("MinTotal: "+str(minTotal)+" was found with angle = "+str(float(minTotalAngle))+"\n")

    FreeCAD.Console.PrintMessage("CSV:\n\n")
    FreeCAD.Console.PrintMessage("Angle,Total Length,Top Line, Right Line, Left Line\n")
    for csv in csvString:
        FreeCAD.Console.PrintMessage(csv)