freecad macro autoload and geometry optimization

fem.py

# Purpose: optimize the thicknesses of the shell to minimize the maximum stress
# The geometry of a shell is set by 12 parameters in the spreadsheet ts.
# Here we directly search for the 12 parameters that minimize the maximum stress,
# subject to the constraint that the volume is 8000 mm^3.
#
# TODO 
# - package as a macro/workbench with buttons
# - 1D option?
# - Z88OS which has "real beam and shell elements" but it's not really possible https://github.com/FreeCAD/FreeCAD/issues/8559
# - maximum is not a good objective. Perhaps the goal is to have each segment
#   have the same (maximum) stress. 
# - generate the sketch and thickness constraints from the script
# - refine the number of parameters. IE. have one segment, optimize the thickness, then
#   make two segments, etc.
# - function evaluations are slow, so maybe a gaussian process model can be worth it
# 
## FreeCAD shell.FCStd
## then in the console
## import femloader
## import fem
## fem.cobyla()
import FreeCAD as App
import FreeCADGui as Gui
import MeshPart
import math
from scipy.optimize import root_scalar, minimize

doc = App.ActiveDocument
sh = App.ActiveDocument.getObject("Spreadsheet")
sk = App.ActiveDocument.getObject("Sketch")
rv = App.ActiveDocument.getObject("Revolution")
bo = App.ActiveDocument.getObject("Body")
an = App.ActiveDocument.getObject("Analysis")


tmin = 0.2

## get the float value of Ai (or ValueError)
def shi(i):
    return float(str(sh.get("A"+str(i))).removesuffix("mm"))

## the volume of the shell is the sum of the volumes of the 12 revolutions of
## rectangles. The top and bottom rectangles are always axis-aligned, so they
## cylinders.
## This might not be right where the segments join. 
def volume():
    v = 0
    # 16 mm high cylinder
    r = 6
    v += 2*math.pi*r*shi(7)*16

    # 12 mm high cylinder
    r = 30
    v += 2*math.pi*r*shi(1)*12


    C2 = float(sh.C2) # z distance covered by every segment

    # a rectangle with sides Ai and vmax where the bottom right corner is at (r, z=0)
    # integrate 2pi r du dv, where u goes over the thickness of the rectangle (0 to Ai)
    # and v goes over the width of the rectangle (0 to vmax)
    # r = r0  + u C2/vmax - Bi v/vmax
    for i in range(5):
        Ai = shi(i+2) # A2, A3, A4, A5, A6
        Bi = shi(i+8) # B1, B2, B3, B4, B5
        vmax = math.sqrt(C2**2 + Bi**2)
        v += 2*math.pi* Ai * ((2*r - Bi)* vmax + Ai * C2)
        r -= Bi
    # missing the final segment

    return v

print(f"volume = {volume()}")

def go(y = []):
  sety(y)
  doc.getObject("FEMMeshNetgen").recompute(True)
  an.recompute(True)
  Gui.Selection.clearSelection()
  Gui.Selection.addSelection('shell', "Analysis")
  Gui.activateWorkbench("FemWorkbench")
  Gui.Selection.clearSelection()
  Gui.Selection.addSelection('shell','SolverCcxTools')
  Gui.runCommand('FEM_SolverRun',0)
  doc.recompute()
  obj = doc.getObject("CCX_Results")
  return max( h - l for h,l in zip(obj.PrincipalMax, obj.PrincipalMin))

def gety():
    y = list()
    i = 0
    # get A1, A2 ... up to the point where a ValueError occurs
    # and put them into y
    while True:
        i += 1
        try:
            y.append(shi(i))
        except ValueError:
            return y

## sety(1) sets all thicknesses to 1
## sety([1]) sets the first thickness to 1
def sety(y):
  if isinstance(y, list):
    for i in range(len(y)):
      sh.set("A"+str(i+1), str(y[i]) + "mm")
  else:
      i = 0
      while True:
          i +=1
          try:
              shi(i)
              sh.set("A"+str(i), str(y) + "mm")
          except ValueError:
              break
  sh.recompute()

def resety():
    for i in range(7):
        sh.set("A"+str(i+1), "1mm")
    for i in range(4):
        sh.set("A"+str(i+8), "4mm")
    sh.set("A12", "20mm")
    sh.recompute()

## multiply all thicknesses by the constant to get the right total volume
def normalizeVolume(target):
    ys = gety()[:7]
    def obj(h):
        sety([ max(tmin, y * math.exp(h)) for y in ys ])
        print(f"recomputing with {h}")
        return volume() - target
    root_scalar(obj, rtol=1e-2, xtol=1e-2, maxiter=10, method="toms748", bracket=[-1, 2])

def bcxt(y):
    return 24 - sum(y[8:12])

def projbcxt(y):
    err = bcxt(y)
    if err < 2:
        k = (24-2) / sum(y[8:12])
        y[8:12] = [ e*k for e in y[8:12]]
        sety(list(y))


# try:
#     normalizeVolume(8000)
# except Exception as e:
#     print(f"failed {e}")


# now I can call cobyla 
nobj = 0
def obj(y, nobj):
    print(f"obj called {nobj}th time with y = {y}")
    sety(list(y))
    projbcxt(y)
    normalizeVolume(8000)
    nobj += 1
    return go()

def maxbs(y):
    return (30-6 - 2) - sum(y[8:12]) 

def maxvol(y):
    sety(list(y))
    return 8000 - volume()

bs = [(0.2, 10) for _ in range(len(gety())-1)]
bs.append((5, 50))
def cobyla():
  return minimize(obj, tuple(gety()), args = (nobj,), method="COBYLA", tol=1e-2, options={"maxiter": 200, "rhobeg": 4, "catol": 1e-2, "disp": True},
        constraints = [ { "type" : "ineq", "fun": maxbs},
                    {"type": "ineq", "fun": maxvol } ],
        bounds= bs)

femloader.py

import sys
import os.path
sys.path += ["/usr/lib/python3/dist-packages/"] # why is this necessary for accessing python3-pyinotify?
import pyinotify
import fem
from importlib import reload


def reload_fem(event):
    try:
        reload(fem)
    except Exception as e:
        print(f"failed to reload fem: {e}")
    print(f"reloaded fem from {event.pathname}")

def reloader():
    # check with pynotify
    wm = pyinotify.WatchManager()
    femfile = os.path.join(os.path.dirname(__file__) , "fem.py")
    wm.add_watch(femfile, pyinotify.IN_MODIFY)
    notifier = pyinotify.ThreadedNotifier(wm, reload_fem)
    notifier.start()
    return notifier

notifier = reloader()