Enhancements to Shaft Design Wizard, e.g. display of stresses for three axes and bending curve for shaft

src/Mod/PartDesign/WizardShaft/SegmentFunction.py

@@ -21,6 +21,7 @@
 # ******************************************************************************/
 
 import FreeCAD # just for debug printing to console...
+import numpy as np
 
 class SegmentFunctionSegment:
     "One segment of a segment function"
@@ -40,6 +41,10 @@ class SegmentFunctionSegment:
         #FIXME: 1E-9 is arbitrary here. But since units are in meters, 1E-9 is a nanometer...
         return abs(self.start - xval) < 1E-9
         
+    def isZero(self):
+        #FIXME: 1E-9 is arbitrary here. But since units are in meters, 1E-9 is a nanometer...
+        return abs(self.coefficient) < 1E-5
+
     def value(self, xval):
         if xval < self.start:
             return 0
@@ -51,10 +56,18 @@ class SegmentFunctionSegment:
         
     def negate(self):
         self.coefficient *= -1
+        return self
+        
+    def negated(self):
+        return SegmentFunctionSegment(self.start, self.variable, self.coefficient * -1.0, self.exponent)
+        
+    def __mul__(self,  value):
+        return SegmentFunctionSegment(self.start, self.variable, self.coefficient * value, self.exponent)
 
     def integrate(self):
         self.exponent = self.exponent + 1
         self.coefficient = self.coefficient * 1 / self.exponent
+        return self
 
     def asString(self):
         return "%f * {%s - %f}^%i" % (self.coefficient, self.variable, self.start, self.exponent)
@@ -70,11 +83,38 @@ class SegmentFunction:
         self.segments = []
         self.name = name
         
+    def findSegment(self,  xval):
+        "Find segment valid for the given xval"
+        for s in self.segments:
+            if s.start <= xval:
+                return s
+        return self.segments[len(self.segments)]
+
+    def isZero(self):
+        for s in self.segments:
+            if not s.isZero():
+                return False
+        return True
+
     def negate(self):
         for s in self.segments:
             s.negate()
         return self
         
+    def negated(self):
+        result = SegmentFunction()
+        result.variable = self.variable
+        for s in self.segments:
+            result.segments.append(s.negated())
+        return result
+        
+    def __mul__(self,  value):
+        result = SegmentFunction()
+        result.variable = self.variable
+        for s in self.segments:
+            result.segments.append(s * value)
+        return result
+
     def index(self, xval):
         "Find insert position for start value xval"
         lastStart = 0.0
@@ -91,10 +131,13 @@ class SegmentFunction:
             #if abs(dict[key]) > 1E-9:
             self.segments.append(SegmentFunctionSegment(key, var, dict[key], 0))
             
+    def addSegment(self, st,  coeff,  exp = 0.0):
+        if abs(coeff) > 1E-9:
+            self.segments.insert(self.index(st), SegmentFunctionSegment(st, self.variable, coeff, exp))
+
     def addSegments(self, dict):
         for key in sorted(dict.iterkeys()):
-            if abs(dict[key]) > 1E-9:
+            self.addSegment(key,  dict[key])
-                self.segments.insert(self.index(key), SegmentFunctionSegment(key, self.variable, dict[key], 0))
 
     def setMaxX(self, mx):
         self.maxX = mx
@@ -124,6 +167,7 @@ class SegmentFunction:
         "Integrate all segments with respect to the variable"
         for s in self.segments:
             s.integrate()
+        return self
 
     def integrated(self):
         "Return a copy of self integrated with respect to the variable"
@@ -156,3 +200,206 @@ class SegmentFunction:
                 FreeCAD.Console.PrintMessage(" + ")
         FreeCAD.Console.PrintMessage("\n")
 
+class IntervalFunction:
+    "Function defined in intervals"
+    intervals = [] # vector of tuples (begin, length)
+    values = [] # vector of constant values applicable for this interval
+    
+    def __init__(self):
+        self.intervals = []
+        self.values = []
+    
+    def addInterval(self,  begin,  length,  value):
+        self.intervals.append((begin,  length))
+        self.values.append(value)
+        
+    def value(self,  xval):            
+        for i in range(len(self.intervals)):
+            if xval >= self.intervals[i][0] and xval < self.intervals[i][0] + self.intervals[i][1]:
+                return  self.values[i]
+        return self.values[len(self.values)-1]
+        
+    def lowervalue(self,  xval):
+        return self.value(xval - 1E-8)
+        
+    def index(self,  xval):
+        lastStart = 0.0
+        for i in range(len(self.intervals)):
+            newStart = self.intervals[i][0]
+            if (xval >= lastStart) and (xval < newStart):
+                return i-1
+            lastStart = newStart
+        return len(self.intervals)-1
+
+    def interval(self, xval):
+        "Return interval (begin, length) for this xval"
+        return self.intervals[self.index(xval)]
+        
+    def begin(self, xval):
+        return self.intervals[self.index(xval)][0]
+        
+    def length(self, xval):
+        return self.intervals[self.index(xval)][1]
+
+class StressFunction:
+    "Specialization for segment-wise display of stresses"
+    # The hairy thing about this is that the segments of the segfunc usually do not correspond with the intervals of the intfunc!
+    segfunc = None # The segment function for the force/moment
+    intfunc = None # The divisors, an interval function giving a specific value for each interval
+    name = "sigma"
+    
+    def __init__(self,  f,  i):
+        self.segfunc = f
+        self.intfunc = i
+        
+    def isZero(self):
+        return self.segfunc.isZero()
+        
+    def evaluate(self, maxX, pointsX):
+        # Note: This usually creates a few more points than specified in pointsX
+        offset = (maxX - self.segfunc.segments[0].start) / (pointsX - 1)
+        xvals = set([self.segfunc.segments[0].start + s * offset for s in range(pointsX)])
+        starts = set([self.segfunc.segments[i].start for i in range(len(self.segfunc.segments))])
+        xvals = xvals.union(starts) # Make sure we have a point on each segment start
+        divs = set([self.intfunc.intervals[i][0] for i in range(len(self.intfunc.intervals))])
+        xvals = xvals.union(divs)
+        
+        xresult = []
+        yresult = []
+        for xval in sorted(xvals):            
+            if xval in starts:
+                # create double point at segment border
+                xresult.append(xval)
+                yresult.append(self.segfunc.lowervalue(xval) / self.intfunc.value(xval))
+            if (xval in divs):
+                # create double point at divisor border
+                xresult.append(xval)
+                yresult.append(self.segfunc.value(xval) / self.intfunc.lowervalue(xval))
+            xresult.append(xval)
+            yresult.append(self.segfunc.value(xval) / self.intfunc.value(xval))
+        return (xresult, yresult)
+        
+class TranslationFunction:
+    "Specialization for segment-wise display of translations"
+    tangfunc = None # The segment function for the tangent to the bending line
+    transfunc = None # The segment function for translations of the shaft (the bending line)
+    intfunc = None # The divisors, a vector of tuples (location, divisor)
+    boundaries = {} # The boundary conditions, dictionary of location:[left boundary, right boundary]
+    module = 210000.0
+    name = "w"
+    
+    def __init__(self,  f,  E,  d,  tangents,  translations):
+        if f.isZero():
+            return
+        # Note: Integration has to be segment-wise because the area moment is not constant in different segments. But this only becomes relevant
+        # when boundary conditions are being applied
+        # E I_i w_i'(x) = tangfunc + C_i0
+        self.tangfunc = f.integrated()
+        self.tangfunc.name = "w'"
+        self.tangfunc.output()
+        # E I_i w_i(x) = transfunc + C_i0 x + C_i1
+        self.transfunc = self.tangfunc.integrated() # + C_i0 * x + C_i1 (integration constants for interval number i)
+        self.transfunc.name = "w"
+        self.transfunc.output()
+        self.module = E
+        self.intfunc = d
+        
+        # Solve boundary conditions. There are two types:
+        # External boundary conditions, e.g. a given tangent direction or translation value at a given x-value
+        # Internal boundary conditions, i.e. at the segment borders the tangent direction and translation of the lines must be equal
+        # Note that the relevant boundaries are those of the intfunc (where the area moment of the shaft cross-section changes)
+        # Every interval of the transfunc has two integration constants C_i0 and C_i1 that need to be defined
+        # Matrix of coefficients
+        A = np.zeros(shape = (2 * len(self.intfunc.intervals),  2 * len(self.intfunc.intervals)))
+        # Vector of RHS values
+        b = np.zeros(shape = 2 * len(self.intfunc.intervals))
+        # Current row where coefficients of next equation will be added
+        row = 0
+        
+        # First look at external boundary conditions
+        for bound in tangents:
+            xval = bound[0]
+            tang = bound[1] 
+            i = self.intfunc.index(xval) # index of this segment
+            I_i = self.intfunc.value(xval) # Area moment of this segment              
+            # w_i'(xval) = tang    =>  (tangfunc(xval) + C_i0) / (E * I_i) = tang =>  C_i0  = tang * (E * I_i) - tangfunc(xval)
+            A[row][2 * i] = 1.0
+            b[row] = tang * E * I_i - self.tangfunc.value(xval)
+            row += 1
+        for bound in translations:
+            xval = bound[0]
+            trans = bound[1] 
+            i = self.intfunc.index(xval) # index of this segment
+            I_i = self.intfunc.value(xval) # Area moment of this segment
+            # w_i(xval) = trans    =>  (transfunc(xval) + C_i0 * xval + C_i1) / (E * I_i) = trans =>  xval / (E * I_i) * C_i0 + 1 / (E * I_i) * C_i1 = trans - transfunc(xval) / (E * I_i)
+            A[row][2 * i] = xval / (E * I_i)
+            A[row][2 * i + 1] = 1 / (E * I_i)
+            b[row] = trans - self.transfunc.value(xval) / (E * I_i)
+            row += 1
+        
+        # Now look at internal boundary conditions (n intervals have n-1 common segment boundaries)
+        for i in range(len(self.intfunc.intervals) - 1):
+            x_start = self.intfunc.intervals[i][0]
+            x_end = x_start + self.intfunc.intervals[i][1]
+            I_i = self.intfunc.value(x_start) # Area moment of this segment
+            I_ip1 = self.intfunc.value(x_end)
+            # w_i'(x_end) = w_i+1'(xend)    =>  (tangfunc(x_end) + C_i0) / (E * I_i) = (tangfunc(x_end) * C_i+1,0) / (E * I_i+1)
+            #   => 1 / (E * I_i) C_i0 - 1 / (E * I_i+1) * C_i+1,0 = tangfunc(x_end) / (E * I_i+1) - tangfunc(x_end) / (E * I_i)
+            A[row][2 * i] = 1 / (E * I_i)
+            A[row][2 * (i+1)] = -1 / (E * I_ip1)
+            b[row] = self.tangfunc.value(x_end) / (E * I_ip1) - self.tangfunc.value(x_end) / (E * I_i)
+            row += 1
+            # w_i(x_end) = w_i+1(xend)    =>  (transfunc(x_end) + C_i0 * x_end + C_i1) / (E * I_i) = (transfunc(x_end) * C_i+1,0) * x_end + C_i+1,1) / (E * I_i+1)
+            #   => x_end / (E * I_i) C_i0 + 1 / (E * I_i) C_i1 - x_end / (E * I_i+1) * C_i+1,0 - 1 / (E * I_i+1) * C_i+1,1 = transfunc(x_end) / (E * I_i+1) - transfunc(x_end) / (E * I_i)
+            A[row][2 * i] = x_end / (E * I_i)
+            A[row][2 * i + 1] = 1 / (E * I_i)
+            A[row][2 * (i+1)] = -x_end / (E * I_ip1)
+            A[row][2 * (i+1) + 1] = -1 / (E * I_ip1)
+            b[row] = self.transfunc.value(x_end) / (E * I_ip1) - self.transfunc.value(x_end) / (E * I_i)
+            row += 1
+            
+        #FreeCAD.Console.PrintMessage(A)
+        #FreeCAD.Console.PrintMessage(" * x = ")
+        #FreeCAD.Console.PrintMessage(b)
+        #FreeCAD.Console.PrintMessage("\n")
+        
+        try:
+            self.boundaries = np.linalg.solve(A, b) # A * self.boundaries = b
+        except np.linalg.linalg.LinAlgError,  e:
+            FreeCAD.Console.PrintMessage(e.message)
+            FreeCAD.Console.PrintMessage(". No solution possible.\n")
+            return
+        
+    def isZero(self):
+        if self.transfunc is None:
+            return True
+        return self.transfunc.isZero()
+        
+    def evaluate(self, maxX, pointsX):
+        # Note: This usually creates a few more points than specified in pointsX
+        offset = (maxX - self.transfunc.segments[0].start) / (pointsX - 1)
+        xvals = set([self.transfunc.segments[0].start + s * offset for s in range(pointsX)])
+        starts = set([self.transfunc.segments[i].start for i in range(len(self.transfunc.segments))])
+        xvals = xvals.union(starts) # Make sure we have a point on each segment start
+        divs = set([self.intfunc.intervals[i][0] for i in range(len(self.intfunc.intervals))])
+        xvals = xvals.union(divs)
+        E = self.module
+        
+        xresult = []
+        yresult = []
+        for xval in sorted(xvals):
+            if xval in divs:
+                i = self.intfunc.index(xval)
+                (begin,  length) = self.intfunc.interval(xval)
+                I_i = self.intfunc.value(xval)
+                C_i0 = self.boundaries[2 * i]
+                C_i1 = self.boundaries[2 * i + 1]
+                FreeCAD.Console.PrintMessage("Interval %u: %f to %f, I_i: %f, C_i0: %f, C_i1: %f\n" % (i,  begin, length,  I_i,  C_i0,  C_i1))           
+            
+            xresult.append(xval)
+            # w(xval) = (transfunc(xval) + C_i0 * xval + C_i1) / (E * I_i)
+            value = (self.transfunc.value(xval)  + C_i0 * xval + C_i1) / (E * I_i)
+            yresult.append(value)
+            
+        return (xresult, yresult)
+

src/Mod/PartDesign/WizardShaft/Shaft.py

@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 #/******************************************************************************
 # *   Copyright (c)2012 Jan Rheinlaender <jrheinlaender@users.sourceforge.net> *
 # *                                                                            *
@@ -20,18 +21,18 @@
 # *                                                                            *
 # ******************************************************************************/
 
-import FreeCAD, FreeCADGui # FreeCAD just required for debug printing to the console...
+import FreeCAD, FreeCADGui
-from SegmentFunction import SegmentFunction
+from SegmentFunction import SegmentFunction,  IntervalFunction,  StressFunction,  TranslationFunction
 from ShaftFeature import ShaftFeature
 from ShaftDiagram import Diagram
+import math
 
 class ShaftSegment:
     length = 0.0
     diameter = 0.0
     innerdiameter = 0.0
-    loadType = "None"
+    constraintType = "None"
-    loadSize = 0.0
+    constraint = None
-    loadLocation = 0.0
 
     def __init__(self, l, d, di):
         self.length = l
@@ -40,22 +41,47 @@ class ShaftSegment:
 
 class Shaft:
     "The axis of the shaft is always assumed to correspond to the X-axis"
+    parent = None
+    doc = None
     # List of shaft segments (each segment has a different diameter)
     segments = []
-    # The sketch
+    # The feature
-    sketch = 0
+    feature = 0
-    #featureWindow = None
     # The diagrams
     diagrams = {} # map of function name against Diagram object
     # Calculation of shaft
-    Qy = 0 # force in direction of y axis
+    F = [None,  None,  None] # force in direction of [x,y,z]-axis
-    Qz = 0 # force in direction of z axis
+    M = [None,  None,  None] # bending moment around [x,z,y]-axis
-    Mbz = 0 # bending moment around z axis
+    w = [None,  None,  None] # Shaft translation due to bending
-    Mby = 0 # bending moment around y axis
+    sigmaN = [None,  None,  None]  # normal stress in direction of x-axis, shear stress in direction of [y,z]-axis
-    Mtz = 0 # torsion moment around z axis
+    sigmaB = [None,  None,  None]  # # torque stress around x-axis, maximum bending stress in direction of [y,z]-axis
+    # Names (note Qy corresponds with Mz, and Qz with My)
+    Fstr = ["Nx","Qy","Qz"] # Forces
+    Mstr = ["Mx","Mz","My"] # Moments
+    wstr = ["",  "wy",  "wz"] # Translations
+    sigmaNstr = ["sigmax","sigmay","sigmaz"] # Normal/shear stresses
+    sigmaBstr = ["taut","sigmabz", "sigmaby"] # Torsion/bending stresses
+    # For diagram labelling
+    Qstrings = (("Normal force [x]",  "x", "mm", "N_x",  "N"), 
+                        ("Shear force [y]",  "x",  "mm",  "Q_y",  "N"), 
+                        ("Shear force [z]",  "x",  "mm",  "Q_z",  "N"))
+    Mstrings = (("Torque [x]",  "x", "mm", "M_t",  "Nm"), 
+                        ("Bending moment [z]",  "x",  "mm",  "M_{b,z}",  "Nm"), 
+                        ("Bending moment [y]",  "x",  "mm",  "M_{b,y}",  "Nm"))
+    wstrings = (("",  "",  "",  "",  ""), 
+                        ("Translation [y]",  "x",  "mm",  "w_y",  "mm"), 
+                        ("Translation [z]",  "x",  "mm",  "w_z",  "mm"))
+    sigmaNstrings =  (("Normal stress [x]",  "x", "mm", "\sigma_x",  u"N/mm²"), 
+                        ("Shear stress [y]",  "x",  "mm",  "\sigma_y",  u"N/mm²"), 
+                        ("Shear stress [z]",  "x",  "mm",  "\sigma_z",  u"N/mm²"))
+    sigmaBstrings = (("Torque stress [x]",  "x", "mm", "\tau_t",  u"N/mm²"), 
+                        ("Bending stress [z]",  "x",  "mm",  "\sigma_{b,z}",  u"N/mm²"), 
+                        ("Bending stress [y]",  "x",  "mm",  "\sigma_{b,y}",  u"N/mm²"))
 
-    def __init__(self, doc):
+    def __init__(self, parent):
-        self.sketch = ShaftFeature(doc)
+        self.parent = parent
+        self.doc = parent.doc
+        self.feature = ShaftFeature(self.doc)
 
     def getLengthTo(self, index):
         "Get the total length of all segments up to the given one"
@@ -65,41 +91,82 @@ class Shaft:
         return result
 
     def addSegment(self, l, d, di):
-        #print "Adding segment: ", l, " : ", d
         self.segments.append(ShaftSegment(l,d,di))
-        self.sketch.addSegment(l, d, di)
+        self.feature.addSegment(l, d, di)
-        # We don't call equilibrium() here because the new segment has no loads defined yet
+        # We don't call equilibrium() here because the new segment has no constraints defined yet
+        # Fix face reference of fixed segment if it is the last one
+        for i in range(1,  len(self.segments)):
+            if self.segments[i].constraintType is not "Fixed":
+                continue
+            if i == len(self.segments) - 1:
+                self.segments[index].constraint.References = [( self.feature.feature,  "Face%u" % (2 * (index+1) + 1) )]
+            else:
+                # Remove reference since it is now in the middle of the shaft (which is not allowed)
+                self.segments[index].constraint.References = [(None,  "")]           
 
     def updateSegment(self, index, length = None, diameter = None, innerdiameter = None):
         oldLength = self.segments[index].length
-        #print "Old length of ", index, ": ", oldLength, ", new Length: ", length, " diameter: ", diameter
+            
         if length is not None:
             self.segments[index].length = length
         if diameter is not None:
             self.segments[index].diameter = diameter
         if innerdiameter is not None:
             self.segments[index].innerdiameter = innerdiameter
-        self.sketch.updateSegment(index, oldLength, self.segments[index].length,
+        
-                                  self.segments[index].diameter, self.segments[index].innerdiameter)
+        self.feature.updateSegment(index, oldLength, self.segments[index].length, self.segments[index].diameter, self.segments[index].innerdiameter)
         self.equilibrium()
         self.updateDiagrams()
 
-    def updateLoad(self, index, loadType = None, loadSize = None, loadLocation = None):
+    def updateConstraint(self, index, constraintType):
-        if (loadType is not None):
+        if (constraintType is not None):
-            self.segments[index].loadType = loadType
+            # Did the constraint type change?
-        if (loadSize is not None):
+            if (self.segments[index].constraintType != "None") and (self.segments[index].constraintType != constraintType):
-            self.segments[index].loadSize = loadSize
+                self.doc.removeObject(self.segments[index].constraint.Name)
-        if (loadLocation is not None):
+                self.segments[index].constraint = None
-            if (loadLocation >= 0) and (loadLocation <= self.segments[index].length):
+            
-                self.segments[index].loadLocation = loadLocation
+            self.segments[index].constraintType = constraintType   
-            else:
+            
-                # TODO: Show warning
+            # Create constraint if it does not exist yet or has changed
-                FreeCAD.Console.PrintMessage("Load location must be inside segment\n")
+            if self.segments[index].constraint is None:
+                if (constraintType == "Force"):
+                    # TODO: Create a reference point and put the force onto it
+                    constraint = self.doc.addObject("Fem::ConstraintForce","ShaftConstraintForce") 
+                    constraint.Force = 1000.0
+                    self.segments[index].constraint = constraint
+                elif (constraintType == "Fixed"):
+                    # TODO: Use robust reference as soon as it is available for the face
+                    constraint = self.doc.addObject("Fem::ConstraintFixed","ShaftConstraintFixed")
+                    if index == 0:
+                        constraint.References = [( self.feature.feature,  "Face1")]
+                    elif index == len(self.segments) - 1:
+                        constraint.References = [( self.feature.feature,  "Face%u" % (2 * (index+1) + 1) )]
+                    self.segments[index].constraint = constraint
+                elif (constraintType == "Bearing"):
+                    # TODO: Use robust reference as soon as it is available for the cylindrical face reference
+                    constraint = self.doc.addObject("Fem::ConstraintBearing","ShaftConstraintBearing")                    
+                    constraint.References = [( self.feature.feature,  "Face%u" % (2 * (index+1)) )]         
+                    constraint.AxialFree = True
+                    self.segments[index].constraint = constraint                    
+                elif (constraintType == "Pulley"):
+                    constraint= self.doc.addObject("Fem::ConstraintPulley","ShaftConstraintPulley")
+                    constraint.References = [( self.feature.feature,  "Face%u" % (2 * (index+1)) )]               
+                    self.segments[index].constraint = constraint                    
+                elif (constraintType == "Gear"):
+                    constraint = self.doc.addObject("Fem::ConstraintGear","ShaftConstraintGear")       
+                    constraint.References = [( self.feature.feature,  "Face%u" % (2 * (index+1)) )]               
+                    self.segments[index].constraint = constraint                      
             
-        #self.feature.updateForces() graphical representation of the forces
         self.equilibrium()
         self.updateDiagrams()
         
+    def editConstraint(self,  index):
+        if (self.segments[index].constraint is not None):
+            FreeCADGui.activeDocument().setEdit(self.segments[index].constraint.Name)
+        
+    def getConstraint(self,  index):
+        return self.segments[index].constraint
+
     def updateEdge(self, column, start):
         App.Console.PrintMessage("Not implemented yet - waiting for robust references...")
         return
@@ -133,110 +200,368 @@ class Shaft:
         return
 
     def updateDiagrams(self):                            
-        if (self.Qy == 0) or (self.Mbz == 0):
+        for ax in range(3):
-            return
+            if self.F[ax] is not None:
-        if self.Qy.name in self.diagrams:
+                if self.F[ax].name in self.diagrams:
-            # Update diagram
+                    self.diagrams[self.F[ax].name].update(self.F[ax], self.getLengthTo(len(self.segments)) / 1000.0)
-            self.diagrams[self.Qy.name].update(self.Qy, self.getLengthTo(len(self.segments)) / 1000.0)
+            if self.M[ax] is not None:
+                if self.M[ax].name in self.diagrams:
+                    self.diagrams[self.M[ax].name].update(self.M[ax], self.getLengthTo(len(self.segments)) / 1000.0)
+            if self.w[ax] is not None:
+                if self.w[ax].name in self.diagrams:
+                    self.diagrams[self.w[ax].name].update(self.w[ax], self.getLengthTo(len(self.segments)) / 1000.0)
+            if self.sigmaN[ax] is not None:
+                if self.sigmaN[ax].name in self.diagrams:
+                    self.diagrams[self.sigmaN[ax].name].update(self.sigmaN[ax], self.getLengthTo(len(self.segments)) / 1000.0)
+            if self.sigmaB[ax] is not None:
+                if self.sigmaB[ax].name in self.diagrams:
+                    self.diagrams[self.sigmaB[ax].name].update(self.sigmaB[ax], self.getLengthTo(len(self.segments)) / 1000.0)
+
+    def showDiagram(self,  which):                
+        if which in self.Fstr:
+            ax = self.Fstr.index(which)
+            text = self.Qstrings[ax]
+            if self.F[ax] == None:
+                # No data
+                return
+            if self.F[ax].name in self.diagrams:
+                # Diagram is already open, close it again
+                self.diagrams[self.F[ax].name].close()
+                del (self.diagrams[self.F[ax].name])
+                return
+            self.diagrams[self.F[ax].name] = Diagram()
+            self.diagrams[self.F[ax].name].create(text[0], self.F[ax], self.getLengthTo(len(self.segments)) / 1000.0, text[1], text[2], 1000.0, text[3], text[4], 1.0, 10)
+        elif which in self.Mstr:
+            ax = self.Mstr.index(which)
+            text = self.Mstrings[ax]
+            if self.M[ax] == None:
+                # No data
+                return
+            if self.M[ax].name in self.diagrams:
+                # Diagram is already open, close it again
+                self.diagrams[self.M[ax].name].close()
+                del (self.diagrams[self.M[ax].name])
+                return
+            self.diagrams[self.M[ax].name] = Diagram()
+            self.diagrams[self.M[ax].name].create(text[0], self.M[ax], self.getLengthTo(len(self.segments)) / 1000.0, text[1], text[2], 1000.0, text[3], text[4], 1.0, 20)
+        elif which in self.wstr:
+            ax = self.wstr.index(which)
+            text = self.wstrings[ax]
+            if self.w[ax] == None:
+                # No data
+                return
+            if self.w[ax].name in self.diagrams:
+                # Diagram is already open, close it again
+                self.diagrams[self.w[ax].name].close()
+                del (self.diagrams[self.w[ax].name])
+                return
+            self.diagrams[self.w[ax].name] = Diagram()
+            self.diagrams[self.w[ax].name].create(text[0], self.w[ax], self.getLengthTo(len(self.segments)) / 1000.0, text[1], text[2], 1000.0, text[3], text[4], 1.0, 30)
+        elif which in self.sigmaNstr:
+            ax = self.sigmaNstr.index(which)
+            text = self.sigmaNstrings[ax]
+            if self.sigmaN[ax] == None:
+                # No data
+                return
+            if self.sigmaN[ax].name in self.diagrams:
+                # Diagram is already open, close it again
+                self.diagrams[self.sigmaN[ax].name].close()
+                del (self.diagrams[self.sigmaN[ax].name])
+                return
+            self.diagrams[self.sigmaN[ax].name] = Diagram()
+            self.diagrams[self.sigmaN[ax].name].create(text[0], self.sigmaN[ax], self.getLengthTo(len(self.segments)) / 1000.0, text[1], text[2], 1000.0, text[3], text[4], 1.0, 10)
+        elif which in self.sigmaBstr:
+            ax = self.sigmaBstr.index(which)
+            text = self.sigmaBstrings[ax]
+            if self.sigmaB[ax] == None:
+                # No data
+                return
+            if self.sigmaB[ax].name in self.diagrams:
+                # Diagram is already open, close it again
+                self.diagrams[self.sigmaB[ax].name].close()
+                del (self.diagrams[self.sigmaB[ax].name])
+                return
+            self.diagrams[self.sigmaB[ax].name] = Diagram()
+            self.diagrams[self.sigmaB[ax].name].create(text[0], self.sigmaB[ax], self.getLengthTo(len(self.segments)) / 1000.0, text[1], text[2], 1000.0, text[3], text[4], 1.0, 20)
+    
+    def addTo(self,  dict,  location,  value):
+        if location not in dict:
+            dict[location] = value
         else:
-            # Create diagram
+            dict[location] += value
-            self.diagrams[self.Qy.name] = Diagram()
-            self.diagrams[self.Qy.name].create("Shear force", self.Qy, self.getLengthTo(len(self.segments)) / 1000.0, "x", "mm", 1000.0, "Q_y", "N", 1.0, 10)
-        if self.Mbz.name in self.diagrams:
-            # Update diagram
-            self.diagrams[self.Mbz.name].update(self.Mbz, self.getLengthTo(len(self.segments)) / 1000.0)
-        else:
-            # Create diagram
-            self.diagrams[self.Mbz.name] = Diagram()
-            self.diagrams[self.Mbz.name].create("Bending moment", self.Mbz, self.getLengthTo(len(self.segments)) / 1000.0, "x", "mm", 1000.0, "M_{b,z}", "Nm", 1.0, 10)
     
     def equilibrium(self):
         # Build equilibrium equations
-        forces = {0.0:0.0} # dictionary of (location : outer force)
-        moments = {0.0:0.0} # dictionary of (location : outer moment)
-        variableNames = [""] # names of all variables
-        locations = {} # dictionary of (variableName : location)
-        coefficientsFy = [0] # force equilibrium equation
-        coefficientsMbz = [0] # moment equilibrium equation
-
-        for i in range(len(self.segments)):
-            lType = self.segments[i].loadType
-            load = -1 # -1 means unknown (just for debug printing)
-            location = -1
-
-            if lType == "Fixed":
-                # Fixed segment
-                if i == 0:
-                    location = 0
-                    variableNames.append("Fy%u" % i)
-                    coefficientsFy.append(1)
-                    coefficientsMbz.append(0)
-                    variableNames.append("Mz%u" % i)
-                    coefficientsFy.append(0)
-                    coefficientsMbz.append(1) # Force does not contribute because location is zero
-                elif i == len(self.segments) - 1:
-                    location = self.getLengthTo(len(self.segments)) / 1000
-                    variableNames.append("Fy%u" % i)
-                    coefficientsFy.append(1)
-                    coefficientsMbz.append(location)
-                    variableNames.append("Mz%u" % i)
-                    coefficientsFy.append(0)
-                    coefficientsMbz.append(1)
-                else:
-                    # TODO: Better error message
-                    FreeCAD.Console.PrintMessage("Fixed constraint must be at beginning or end of shaft\n")
-                    return
-
-                locations["Fy%u" % i] = location
-                locations["Mz%u" % i] = location
-            elif lType == "Static":
-                # Static load (currently force only)
-                load = self.segments[i].loadSize
-                location = (self.getLengthTo(i) + self.segments[i].loadLocation)  / 1000 # convert to meters
-                coefficientsFy[0] = coefficientsFy[0] - load
-                forces[location] = load
-                coefficientsMbz[0] = coefficientsMbz[0] - load * location
-                moments[location] = 0
-            #elif lType == "None":
-            #    # No loads on segment
-
-            FreeCAD.Console.PrintMessage("Segment: %u, type: %s, load: %f, location: %f\n" % (i, lType, load, location))
-
-        self.printEquilibrium(variableNames, coefficientsFy)
-        self.printEquilibrium(variableNames, coefficientsMbz)
-
-        # Build matrix and vector for linear algebra solving algorithm
         try:
             import numpy as np
         except ImportError:
             FreeCAD.Console.PrintMessage("numpy is not installed on your system\n")
             raise ImportError("numpy not installed")
-        if (len(coefficientsFy) < 3) or (len(coefficientsMbz) < 3):
-            return
-        A = np.array([coefficientsFy[1:], coefficientsMbz[1:]])
-        b = np.array([coefficientsFy[0], coefficientsMbz[0]])
-        solution = np.linalg.solve(A, b)
                 
-        # Complete dictionary of forces and moments
+        # Initialization of structures. All three axes are handled separately so everything is 3-fold
-        if variableNames[1][0] == "F":
+        # dictionaries of (location : outer force/moment) with reverse sign, which means that the segment functions for the section force and section moment
-            forces[locations[variableNames[1]]] = solution[0]
+        # created from them will have signs as by the convention in 
-        else:
+        # http://www.umwelt-campus.de/ucb/fileadmin/users/90_t.preussler/dokumente/Skripte/TEMECH/TMI/Ebene_Balkenstatik.pdf (page 10)
-            moments[locations[variableNames[1]]] = solution[0]
+        # (see also example on page 19)
+        forces = [{0.0:0.0},  {0.0:0.0},  {0.0:0.0}] 
+        moments = [{0.0:0.0},  {0.0:0.0},  {0.0:0.0}]
+        # Boundary conditions for shaft bending line
+        tangents = [[],  [],  []] # Tangents to shaft bending line
+        translations = [[],  [],  []] # Shaft displacement
+        # Variable names, e.g. Fx, Mz. Because the system must be exactly determined, not more than two independent variables for each 
+        # force/moment per axis are possible (if there are more no solution is calculated)
+        variableNames = [[""],  [""], [""]] 
+        # # dictionary of (variableName : location) giving the x-coordinate at which the force/moment represented by the variable acts on the shaft
+        locations = {}
+        # Coefficients of the equilibrium equations in the form a = b * F1 + c * F2 and d = e * M1 + f * M2
+        # LHS (variables a1, a2, a3, d3) initialized to zero
+        coefficientsF = [[0], [0], [0]]
+        coefficientsM = [[0], [0], [0]]
 
-        if variableNames[2][0] == "F":
+        for i in range(len(self.segments)):
-            forces[locations[variableNames[2]]] = solution[1]
+            cType = self.segments[i].constraintType
-        else:
+            constraint = self.segments[i].constraint
-            moments[locations[variableNames[2]]] = solution[1]
             
-        FreeCAD.Console.PrintMessage(forces)
+            if cType == "Fixed":
-        FreeCAD.Console.PrintMessage(moments)
+                # Fixed segment
-        self.Qy = SegmentFunction("Qy")
+                if i == 0:
-        self.Qy.buildFromDict("x", forces)
+                    # At beginning of shaft
-        self.Qy.output()
+                    location = 0                    
-        self.Mbz = self.Qy.integrated().negate()
+                elif i == len(self.segments) - 1:
-        self.Mbz.addSegments(moments) # takes care of boundary conditions
+                    # At end of shaft
-        self.Mbz.name = "Mbz"
+                    location = self.getLengthTo(len(self.segments)) / 1000.0 # convert to meters                    
-        self.Mbz.output()
+                else:
+                    # TODO: Better error message
+                    FreeCAD.Console.PrintMessage("Fixed constraint must be at beginning or end of shaft\n")
+                    return
+                    
+                for ax in range(3):        
+                    # Create a new reaction force
+                    variableNames[ax].append("%s%u" % (self.Fstr[ax], i))
+                    coefficientsF[ax].append(1)
+                    # Register location of reaction force             
+                    locations["%s%u" % (self.Fstr[ax], i)] = location
+                    # Boundary conditions for the translations
+                    tangents[ax].append((location,  0.0))
+                    translations[ax].append((location,  0.0))
+                coefficientsM[0].append(0) # Reaction force contributes no moment around x axis
+                coefficientsM[1].append(location) # Reaction force contributes a positive moment around z axis
+                coefficientsM[2].append(-location) # Reaction force contributes a negative moment around y axis  
+                
+                for ax in range(3):
+                    # Create a new reaction moment
+                    variableNames[ax].append("%s%u" % (self.Mstr[ax], i))
+                    coefficientsF[ax].append(0)
+                    coefficientsM[ax].append(1)                                   
+                    locations["%s%u" % (self.Mstr[ax], i)] = location
+                    
+            elif cType == "Force":
+                # Static force (currently force on midpoint of segment only)
+                force = constraint.DirectionVector.multiply(constraint.Force)                
+                # TODO: Extract value of the location from geometry
+                location = (self.getLengthTo(i) + self.segments[i].length/2.0)  / 1000.0
+                # The force itself
+                for ax in range(3): 
+                    if abs(force[ax]) > 0.0:
+                        coefficientsF[ax][0] = coefficientsF[ax][0] - force[ax] # neg. because this coefficient is on the LHS of the equilibrium equation
+                        self.addTo(forces[ax],  location,  -force[ax]) # neg. to fulfill the convention mentioned above
+                # Moments created by the force (by definition no moment is created by the force in x-direction)                
+                if abs(force[1]) > 0.0:
+                    coefficientsM[1][0] = coefficientsM[1][0] - force[1] * location # moment around z-axis
+                    self.addTo(moments[1],  location,  0)
+                if abs(force[2]) > 0.0:
+                    coefficientsM[2][0] = coefficientsM[2][0] + force[2] * location # moment around y-axis
+                    self.addTo(moments[2],  location,  0) # No outer moment acts here!
+                    
+            elif cType == "Bearing":
+                 location = constraint.BasePoint.x / 1000.0 # TODO: This assumes that the shaft feature starts with the first segment at (0,0,0)  and its axis corresponds to the x-axis
+                 # Bearing reaction forces. TODO: the bearing is assumed to not induce any reaction moments
+                 start = (0 if constraint.AxialFree == False else 1)
+                 for ax in range(start, 3):
+                    variableNames[ax].append("%s%u" % (self.Fstr[ax], i))
+                    coefficientsF[ax].append(1)
+                    locations["%s%u" % (self.Fstr[ax], i)] = location
+                    # Boundary condition
+                    translations[ax].append((location,  0.0))
+                 if constraint.AxialFree == False:
+                    coefficientsM[0].append(0)             # Reaction force contributes no moment around x axis
+                 coefficientsM[1].append(location) # Reaction force contributes a positive moment around z axis
+                 coefficientsM[2].append(-location) # Reaction force contributes a negative moment around y axis                 
+                
+            elif cType == "Gear":
+                force = constraint.DirectionVector.multiply(constraint.Force) 
+                location = constraint.BasePoint.x / 1000.0
+                lever = [0,  constraint.Diameter/2.0/1000.0 * math.sin(constraint.ForceAngle / 180.0 * math.pi),  
+                                    constraint.Diameter/2.0 /1000.0* math.cos(constraint.ForceAngle / 180.0 * math.pi)]
+                
+                # Effect of the gear force
+                for ax in range(3):
+                    if abs(force[ax]) > 0.0:
+                        # Effect of the force
+                        coefficientsF[ax][0] = coefficientsF[ax][0] - force[ax]
+                        self.addTo(forces[ax],  location,  -force[ax])
+                # Moments created by the force (by definition no moment is created by the force in x-direction)                
+                if abs(force[1]) > 0.0:
+                    coefficientsM[1][0] = coefficientsM[1][0] - force[1] * location # moment around z-axis
+                    self.addTo(moments[1],  location,  0)
+                if abs(force[2]) > 0.0:
+                    coefficientsM[2][0] = coefficientsM[2][0] + force[2] * location # moment around y-axis
+                    self.addTo(moments[2],  location,  0) # No outer moment acts here!
+                
+                # Moments created by the force and lever
+                if abs(force[0]) > 0.0:
+                    momenty = force[0] * lever[2]
+                    momentz = force[0] * lever[1]
+                    coefficientsM[1][0] = coefficientsM[1][0] + momentz # moment around z-axis
+                    self.addTo(moments[1],  location,  momentz)
+                    coefficientsM[2][0] = coefficientsM[2][0] - momenty # moment around y-axis
+                    self.addTo(moments[2],  location,  -momenty)         
+                if abs(force[1]) > 0.0:
+                    moment = force[1] * lever[2]
+                    coefficientsM[0][0] = coefficientsM[0][0] + moment
+                    self.addTo(moments[0],  location,  moment)
+                if abs(force[2]) > 0.0:
+                    moment = force[2] * lever[1]
+                    coefficientsM[0][0] = coefficientsM[0][0] - moment
+                    self.addTo(moments[0],  location,  -moment)
+            elif cType == "Pulley":
+                forceAngle1 = (constraint.ForceAngle + constraint.BeltAngle + 90.0) / 180.0 * math.pi
+                forceAngle2 = (constraint.ForceAngle - constraint.BeltAngle + 90.0) / 180.0 * math.pi
+                #FreeCAD.Console.PrintMessage("BeltForce1: %f, BeltForce2: %f\n" % (constraint.BeltForce1,  constraint.BeltForce2))
+                #FreeCAD.Console.PrintMessage("Angle1: %f, Angle2: %f\n" % (forceAngle1,  forceAngle2))
+                force = [0,  -constraint.BeltForce1 * math.sin(forceAngle1) - constraint.BeltForce2 * math.sin(forceAngle2), 
+                                   constraint.BeltForce1 * math.cos(forceAngle1) + constraint.BeltForce2 * math.cos(forceAngle2)]
+                location = constraint.BasePoint.x / 1000.0                
+                
+                # Effect of the pulley forces
+                for ax in range(3):
+                    if abs(force[ax]) > 0.0:
+                        # Effect of the force
+                        coefficientsF[ax][0] = coefficientsF[ax][0] - force[ax]
+                        self.addTo(forces[ax],  location,  -force[ax])
+                # Moments created by the force (by definition no moment is created by the force in x-direction)                
+                if abs(force[1] ) > 0.0:
+                    coefficientsM[1][0] = coefficientsM[1][0] - force[1] * location # moment around z-axis
+                    self.addTo(moments[1],  location,  0)
+                if abs(force[2]) > 0.0:
+                    coefficientsM[2][0] = coefficientsM[2][0] + force[2] * location # moment around y-axis
+                    self.addTo(moments[2],  location,  0) # No outer moment acts here!
+                
+                # Torque
+                moment = constraint.Force * (1 if constraint.IsDriven is True else -1) 
+                coefficientsM[0][0] = coefficientsM[0][0] + moment
+                self.addTo(moments[0],  location,  moment)
+
+        areas = [None,  None,  None]
+        areamoments = [None,  None,  None]
+        torquemoments = [None,  None,  None]
+        
+        for ax in range(3):
+            FreeCAD.Console.PrintMessage("Axis: %u\n" %  ax) 
+            self.printEquilibrium(variableNames[ax], coefficientsF[ax])
+            self.printEquilibrium(variableNames[ax], coefficientsM[ax])
+
+            if  len(coefficientsF[ax]) <= 1:
+                # Note: coefficientsF and coefficientsM always have the same length
+                FreeCAD.Console.PrintMessage("Matrix is singular, no solution possible\n")
+                self.parent.updateButtons(ax,  False)
+                continue
+                
+            # Handle special cases. Note that the code above should ensure that coefficientsF and coefficientsM always have same length
+            solution = [None,  None]
+            if len(coefficientsF[ax]) == 2:
+                if coefficientsF[ax][1] != 0.0 and coefficientsF[ax][0] != 0.0:
+                    solution[0] = coefficientsF[ax][0] / coefficientsF[ax][1]
+                if coefficientsM[ax][1] != 0.0 and coefficientsM[ax][0] != 0.0:
+                    solution[1] = coefficientsM[ax][0] / coefficientsM[ax][1]
+                    if abs(solution[0] - solution[1]) < 1E9:
+                        FreeCAD.Console.PrintMessage("System is statically undetermined. No solution possible.\n")
+                        self.parent.updateButtons(ax,  False)
+                        continue
+            else:            
+                # Build matrix and vector for linear algebra solving algorithm           
+               # TODO: This could easily be done manually... there are only 2 variables and 6 coefficients
+                A = np.array([coefficientsF[ax][1:], coefficientsM[ax][1:]])                           
+                b = np.array([coefficientsF[ax][0], coefficientsM[ax][0]])
+                try:
+                    solution = np.linalg.solve(A, b) # A * solution = b
+                except np.linalg.linalg.LinAlgError,  e:
+                    FreeCAD.Console.PrintMessage(e.message)
+                    FreeCAD.Console.PrintMessage(". No solution possible.\n")
+                    self.parent.updateButtons(ax,  False)
+                    continue
+
+            # Complete dictionary of forces and moments with the two reaction forces that were calculated
+            for i in range(2):
+                if solution[i] is None:
+                    continue
+                FreeCAD.Console.PrintMessage("Reaction force/moment: %s = %f\n" % (variableNames[ax][i+1],  solution[i]))
+                if variableNames[ax][i+1][0] == "M":
+                    moments[ax][locations[variableNames[ax][i+1]]] = -solution[i]
+                else:
+                    forces[ax][locations[variableNames[ax][i+1]]] = -solution[i]
+
+            FreeCAD.Console.PrintMessage(forces[ax])
+            FreeCAD.Console.PrintMessage("\n")
+            FreeCAD.Console.PrintMessage(moments[ax])
+            FreeCAD.Console.PrintMessage("\n")
+            
+            # Forces
+            self.F[ax] = SegmentFunction(self.Fstr[ax])
+            self.F[ax].buildFromDict("x", forces[ax])
+            self.parent.updateButton(1,  ax,  not self.F[ax].isZero())
+            self.F[ax].output()        
+            # Moments
+            if ax == 0:
+                self.M[0] = SegmentFunction(self.Mstr[0])
+                self.M[0].buildFromDict("x",  moments[0])
+            elif ax == 1:
+                self.M[1] = self.F[1].integrated().negate()
+                self.M[1].name = self.Mstr[1]
+                self.M[1].addSegments(moments[1]) # takes care of boundary conditions       
+            elif ax == 2:
+                self.M[2] = self.F[2].integrated()        
+                self.M[2].name = self.Mstr[2]
+                self.M[2].addSegments(moments[2]) # takes care of boundary conditions       
+            self.parent.updateButton(2,  ax,  not self.M[ax].isZero())
+            self.M[ax].output()
+            
+            # Areas and area moments            
+            location = 0.0
+            areas[ax] = IntervalFunction()
+            areamoments[ax] = IntervalFunction()
+            torquemoments[ax] = IntervalFunction()
+            
+            for i in range(len(self.segments)):
+                od = self.segments[i].diameter
+                id = self.segments[i].innerdiameter
+                length = self.segments[i].length/1000.0
+                areas[ax].addInterval(location, length,  math.pi/4.0 * (math.pow(self.segments[i].diameter,  2.0) - math.pow(self.segments[i].innerdiameter,  2.0)))
+                areamoment = math.pi/32.0 * (math.pow(self.segments[i].diameter,  4.0) - math.pow(self.segments[i].innerdiameter,  4.0)) / self.segments[i].diameter
+                areamoments[ax].addInterval(location,  length,  areamoment)
+                torquemoments[ax].addInterval(location,  length,  2 * areamoment)
+                location += length
+            
+            # Bending line
+            if ax > 0:
+                if len(tangents[ax])+ len(translations[ax]) == 2:
+                    # TODO: Get Young's module from material type instead of using 210000 N/mm²
+                    self.w[ax] = TranslationFunction(self.M[ax].negated() * 1000.0,  210000,  areamoments[ax], tangents[ax], translations[ax])
+                    self.w[ax].name= self.wstr[ax]
+                    self.parent.updateButton(3,  ax,  not self.w[ax].isZero())
+                else:
+                    self.parent.updateButton(3,  ax,  False)                
+            
+            # Normal/shear stresses and torque/bending stresses
+            self.sigmaN[ax] = StressFunction(self.F[ax],  areas[ax])
+            self.sigmaN[ax].name = self.sigmaNstr[ax]
+            self.parent.updateButton(4,  ax,  not self.sigmaN[ax].isZero())
+            if ax == 0:
+                self.sigmaB[ax] = StressFunction(self.M[ax] * 1000.0,  torquemoments[ax])
+            else:
+                self.sigmaB[ax] = StressFunction(self.M[ax] * 1000.0,  areamoments[ax])
+            self.sigmaB[ax].name = self.sigmaBstr[ax]
+            self.parent.updateButton(5,  ax,  not self.sigmaB[ax].isZero())
 
     def printEquilibrium(self, var, coeff):
         # Auxiliary method for debugging purposes        

src/Mod/PartDesign/WizardShaft/ShaftDiagram.py

@@ -38,6 +38,7 @@ class Diagram:
     ypoints = []
     # Plot object
     thePlot = None
+    win = None
 
     def create(self, title, function, xlength, xname, xunit, xscale, yname, yunit, yscale, numxpoints):
         # Initialize
@@ -54,7 +55,7 @@ class Diagram:
         self.numxpoints = numxpoints
 
         # Create a plot window
-        win = Plot.figure(title)
+        self.win = Plot.figure(title)
         # Get the plot object from the window
         self.thePlot = Plot.getPlot()
         # Format the plot object
@@ -96,3 +97,7 @@ class Diagram:
         axes.set_xlim(right = max(self.xpoints) * 1.05)
         axes.set_ylim(min(self.ypoints) * 1.05, max(self.ypoints) * 1.05)
         self.thePlot.update()
+        
+    def close(self):
+        # Close the associated mdiSubWindow
+        self.win.parent().close()

src/Mod/PartDesign/WizardShaft/WizardShaft.py

@@ -21,7 +21,6 @@
 # ******************************************************************************/
 
 import FreeCAD, FreeCADGui
-#import os
 from PyQt4 import QtCore, QtGui
 from WizardShaftTable import WizardShaftTable
 from Shaft import Shaft
@@ -38,6 +37,8 @@ class TaskWizardShaft:
     shaft = 0
     # Feature
     featureWindow = 0
+    # Buttons
+    buttons = [[None,  None,  None],  [None,  None,  None],  [None,  None,  None],  [None,  None,  None],  [None,  None,  None],  [None,  None,  None]]
 
     def __init__(self, doc):
         mw = QtGui.qApp.activeWindow()
@@ -55,13 +56,95 @@ class TaskWizardShaft:
         else:
             featureWindow = cw.activeSubWindow()
         
-        # Create Shaft object
+        # Buttons for diagram display
-        self.shaft = Shaft(self.doc)
+        buttons = QtGui.QGridLayout() 
+        bnames = [["All [x]", "All [y]", "All [z]" ], 
+                           ["N [x]", "Q [y]", "Q [z]"], 
+                           ["Mt [x]",  "Mb [z]", "Mb [y]"], 
+                           ["",  "w [y]",  "w [z]"], 
+                           ["sigma [x]", "sigma [y]", "sigma [z]"], 
+                           ["tau [x]",  "sigmab [z]", "sigmab [y]"]]
+        slots = [[self.slotAllx,  self.slotAlly,  self.slotAllz], 
+                      [self.slotFx,  self.slotQy,  self.slotQz], 
+                      [self.slotMx,  self.slotMz,  self.slotMy], 
+                      [self.slotNone,  self.slotWy, self.slotWz], 
+                      [self.slotSigmax,  self.slotSigmay,  self.slotSigmaz], 
+                      [self.slotTaut,  self.slotSigmabz,  self.slotSigmaby]]
+        for col in range(3):
+            for row in range(6):
+                button = QtGui.QPushButton(bnames[row][col])
+                buttons.addWidget(button,  row,  col)
+                self.buttons[row][col] = button
+                button.clicked.connect(slots[row][col])
                 
-        # Assign a table widget to the dock window
+        # Create Shaft object
+        self.shaft = Shaft(self)
+        # Create table widget
+        self.form = QtGui.QWidget()
         self.table = WizardShaftTable(self, self.shaft)        
-        self.form = self.table.widget
+        
-        self.form.setWindowTitle("Shaft wizard")
+        # The top layout will contain the Shaft Wizard layout plus the elements of the FEM constraints dialog
+        layout = QtGui.QVBoxLayout()
+        layout.setObjectName("ShaftWizard") # Do not change or translate: Required to detect whether Shaft Wizard is running in FemGui::ViewProviderFemConstraintXXX
+        sublayout = QtGui.QVBoxLayout()
+        sublayout.setObjectName("ShaftWizardLayout") # Do not change or translate        
+        sublayout.addWidget(self.table.widget)
+        sublayout.addLayout(buttons)
+        layout.addLayout(sublayout)
+        self.form.setLayout(layout)
+        
+        # Switch to feature window
+        mdi=QtGui.qApp.activeWindow().findChild(QtGui.QMdiArea)
+        cw.setActiveSubWindow(featureWindow)
+        
+    def slotAllx(self):
+        self.shaft.showDiagram("Allx")
+    def slotAlly(self):
+        self.shaft.showDiagram("Ally")
+    def slotAllz(self):
+        self.shaft.showDiagram("Allz")
+        
+    def slotFx(self):
+        self.shaft.showDiagram("Nx")        
+    def slotQy(self):
+        self.shaft.showDiagram("Qy")        
+    def slotQz(self):
+        self.shaft.showDiagram("Qz")
+    
+    def slotMx(self):
+        self.shaft.showDiagram("Mx")
+    def slotMz(self):
+        self.shaft.showDiagram("Mz")        
+    def slotMy(self):
+        self.shaft.showDiagram("My")
+        
+    def slotNone(self):
+        pass
+    def slotWy(self):
+        self.shaft.showDiagram("wy")
+    def slotWz(self):
+        self.shaft.showDiagram("wz")
+        
+    def slotSigmax(self):
+        self.shaft.showDiagram("sigmax")
+    def slotSigmay(self):
+        self.shaft.showDiagram("sigmay")        
+    def slotSigmaz(self):
+        self.shaft.showDiagram("sigmaz")
+        
+    def slotTaut(self):
+        self.shaft.showDiagram("taut")
+    def slotSigmabz(self):
+        self.shaft.showDiagram("sigmabz")        
+    def slotSigmaby(self):
+        self.shaft.showDiagram("sigmaby")
+        
+    def updateButton(self,  row,  col,  flag):
+        self.buttons[row][col].setEnabled(flag)
+        
+    def updateButtons(self,  col,  flag):
+        for row in range(len(self.buttons)):
+            self.updateButton(row,  col,  flag)
         
     def getStandardButtons(self):
         return int(QtGui.QDialogButtonBox.Ok)
@@ -75,9 +158,19 @@ class TaskWizardShaft:
             del self.form
         return True
 
+# Work-around to allow a callback
+# Problem: From the FemConstraint ViewProvider, we need to tell the Shaft instance that the user finished editing the constraint
+# We can find the Shaft Wizard dialog object from C++, but there is not way to reach the Shaft instance
+# Also it seems to be impossible to access the active dialog from Python, so Gui::Command::runCommand() is not an option either
+# Note: Another way would be to create a hidden widget in the Shaft Wizard dialog and write some data to it, triggering a slot 
+# in the python code
+WizardShaftDlg = None
+
 class WizardShaftGui:    
     def Activated(self):
-        FreeCADGui.Control.showDialog(TaskWizardShaft(FreeCAD.ActiveDocument))
+        global WizardShaftDlg
+        WizardShaftDlg = TaskWizardShaft(FreeCAD.ActiveDocument)
+        FreeCADGui.Control.showDialog(WizardShaftDlg)
         
     def GetResources(self):
         IconPath = FreeCAD.ConfigGet("AppHomePath") + "Mod/PartDesign/WizardShaft/WizardShaft.svg"
@@ -88,7 +181,28 @@ class WizardShaftGui:
     def IsActive(self):
         return FreeCAD.ActiveDocument != None
         
+    def __del__(self):
+        global WizardShaftDlg
+        WizardShaftDlg = None
+        
+class WizardShaftGuiCallback:    
+    def Activated(self):
+        global WizardShaftDlg
+        if WizardShaftDlg != None and WizardShaftDlg.table != None:
+            WizardShaftDlg.table.finishEditConstraint()
+                
+    def isActive(self):
+        global WizardShaftDlg
+        return (WizardShaftDlg is not None)
+        
+    def GetResources(self):
+        IconPath = FreeCAD.ConfigGet("AppHomePath") + "Mod/PartDesign/WizardShaft/WizardShaft.svg"
+        MenuText = 'Shaft design wizard...'
+        ToolTip  = 'Start the shaft design wizard'
+        return {'Pixmap' : IconPath, 'MenuText': MenuText, 'ToolTip': ToolTip}
+
 FreeCADGui.addCommand('PartDesign_WizardShaft', WizardShaftGui())
+FreeCADGui.addCommand('PartDesign_WizardShaftCallBack', WizardShaftGuiCallback())
 
 #Note: Start wizard in Python Console with
 # Gui.runCommand('PartDesign_WizardShaft')

src/Mod/PartDesign/WizardShaft/WizardShaftTable.py

@@ -31,21 +31,18 @@ class WizardShaftTable:
         "Length"        : 0,
         "Diameter"      : 1,
         "InnerDiameter" : 2,
-        "LoadType"      : 3,
+        "ConstraintType"      : 3,
-        "LoadSize"      : 4,
+        "StartEdgeType" : 4,
-        "LoadLocation"  : 5,
+        "StartEdgeSize" : 5,
-        "StartEdgeType" : 6,
+        "EndEdgeType"   : 6,
-        "StartEdgeSize" : 7,
+        "EndEdgeSize"   : 7
-        "EndEdgeType"   : 8,
-        "EndEdgeSize"   : 9
     }
     rowDictReverse = {}
-    headers = ["Length [mm]",
+    headers = [
+                "Length [mm]",
                "Diameter [mm]",
                "Inner diameter [mm]",
-               "Load type",
+               "Constraint type",
-               "Load [N]",
-               "Location [mm]",
                "Start edge type",
                "Start edge size",
                "End edge type",
@@ -54,6 +51,8 @@ class WizardShaftTable:
     widget = 0
     wizard = 0
     shaft = 0
+    editedRow = None
+    editedColumn = None
 
     def __init__(self, w, s):
         for key in self.rowDict.iterkeys():
@@ -63,8 +62,9 @@ class WizardShaftTable:
         self.shaft = s
         # Create table widget
         self.widget = QtGui.QTableWidget(len(self.rowDict), 0)     
+        self.widget.setObjectName("ShaftWizardTable") # Do not change or translate: Used in ViewProviderFemConstraintXXX
+        self.widget.setWindowTitle("Shaft wizard")
         self.widget.resize(QtCore.QSize(300,200))
-        #self.widget.setFocusPolicy(QtCore.Qt.StrongFocus)
 
         # Label rows and columns
         self.widget.setVerticalHeaderLabels(self.headers)
@@ -82,14 +82,12 @@ class WizardShaftTable:
         self.addColumn()
         self.setLength(0, 40.0)
         self.setDiameter(0, 50.0)
-        self.setLoadType(0, "Static")
+        self.setConstraintType(0, "Bearing")
-        self.setLoadSize(0, 1000.0)
-        self.setLoadLocation(0, 25.0)
         # Section 2
         self.addColumn()
         self.setLength(1, 80.0)
         self.setDiameter(1, 60.0)
-        self.setLoadType(1, "Fixed")
+        self.setConstraintType(1, "Force")
 
     def slotInsertColumn(self, point):
         # FIXME: Allow inserting columns, not just adding at the end
@@ -144,32 +142,21 @@ class WizardShaftTable:
         widget.setValue(innerdiameter)
         widget.valueChanged.connect(self.slotValueChanged)
         widget.editingFinished.connect(self.slotEditingFinished)
-        # Load type
+        # Constraint type
         widget = QtGui.QComboBox(self.widget)
         widget.insertItem(0, "None")
         widget.insertItem(1, "Fixed")
-        widget.insertItem(2, "Static")
+        widget.insertItem(2, "Force")
         widget.insertItem(3, "Bearing")
-        widget.insertItem(4, "Pulley")
+        widget.insertItem(4, "Gear")
-        self.widget.setCellWidget(self.rowDict["LoadType"], index, widget)
+        widget.insertItem(5, "Pulley")
+        action = QtGui.QAction("Edit constraint", widget)
+        action.triggered.connect(self.slotEditConstraint)
+        widget.addAction(action)
+        widget.setContextMenuPolicy(QtCore.Qt.ActionsContextMenu)
+        self.widget.setCellWidget(self.rowDict["ConstraintType"], index, widget)
         widget.setCurrentIndex(0)
-        self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotLoadType)
+        self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotConstraintType)      
-        # Load size
-        widget = QtGui.QDoubleSpinBox(self.widget)
-        widget.setMinimum(-1E9)
-        widget.setMaximum(1E9)
-        self.widget.setCellWidget(self.rowDict["LoadSize"], index, widget)
-        widget.setValue(0)
-        widget.valueChanged.connect(self.slotValueChanged)
-        widget.editingFinished.connect(self.slotEditingFinished)
-        # Load location
-        widget = QtGui.QDoubleSpinBox(self.widget)
-        widget.setMinimum(0)
-        widget.setMaximum(1E9)
-        self.widget.setCellWidget(self.rowDict["LoadLocation"], index, widget)
-        widget.setValue(0)
-        widget.valueChanged.connect(self.slotValueChanged)
-        widget.editingFinished.connect(self.slotEditingFinished)
         # Start edge type
         widget = QtGui.QComboBox(self.widget)
         widget.insertItem(0, "None",)
@@ -177,7 +164,7 @@ class WizardShaftTable:
         widget.insertItem(2, "Fillet")
         self.widget.setCellWidget(self.rowDict["StartEdgeType"],index, widget)
         widget.setCurrentIndex(0)
-        self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotLoadType)
+        #self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotLoadType)
         # Start edge size
         widget = QtGui.QDoubleSpinBox(self.widget)
         widget.setMinimum(0)
@@ -193,7 +180,7 @@ class WizardShaftTable:
         widget.insertItem(2, "Fillet")
         self.widget.setCellWidget(self.rowDict["EndEdgeType"],index, widget)
         widget.setCurrentIndex(0)
-        self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotLoadType)
+        #self.widget.connect(widget, QtCore.SIGNAL("currentIndexChanged(const QString&)"), self.slotLoadType)
         # End edge size
         widget = QtGui.QDoubleSpinBox(self.widget)
         widget.setMinimum(0)
@@ -208,6 +195,8 @@ class WizardShaftTable:
         self.editedValue = value
 
     def slotEditingFinished(self):
+        if self.editedRow == None:
+            return
         rowName = self.rowDictReverse[self.editedRow]
         if rowName is None:
             return
@@ -217,12 +206,8 @@ class WizardShaftTable:
             self.shaft.updateSegment(self.editedColumn, diameter = self.getDoubleValue(rowName, self.editedColumn))
         elif rowName == "InnerDiameter":
             self.shaft.updateSegment(self.editedColumn, innerdiameter = self.getDoubleValue(rowName, self.editedColumn))
-        elif rowName == "LoadType":
+        elif rowName == "Constraintype":
-            self.shaft.updateLoad(self.editedColumn, loadType = self.getListValue(rowName, self.editedColumn))
+            self.shaft.updateConstraint(self.editedColumn, self.getListValue(rowName, self.editedColumn))
-        elif rowName == "LoadSize":
-            self.shaft.updateLoad(self.editedColumn, loadSize = self.getDoubleValue(rowName, self.editedColumn))
-        elif rowName == "LoadLocation":
-            self.shaft.updateLoad(self.editedColumn, loadLocation = self.getDoubleValue(rowName, self.editedColumn))
         elif rowName == "StartEdgeType":
             pass
         elif rowName == "StartEdgeSize":
@@ -232,6 +217,13 @@ class WizardShaftTable:
         elif rowName == "EndEdgeSize":
             pass
 
+    def slotEditConstraint(self):        
+        (self.editedRow, self.editedColumn) = self.getFocusedCell() # Because finishEditConstraint() will trigger slotEditingFinished() which requires this information
+        self.shaft.editConstraint(self.editedColumn)
+        
+    def finishEditConstraint(self):
+        self.shaft.updateConstraint(self.editedColumn, self.getConstraintType(self.editedColumn))
+        
     def setLength(self, column, l):
         self.setDoubleValue("Length", column, l)
         self.shaft.updateSegment(column, length = l)
@@ -254,32 +246,15 @@ class WizardShaftTable:
         return self.getDoubleValue("InnerDiameter", column)
 
     @QtCore.pyqtSlot('QString')
-    def slotLoadType(self, text):
+    def slotConstraintType(self, text):
-        if text != "Fixed":
+        self.shaft.updateConstraint(self.getFocusedColumn(), text)
-            if (self.getLoadSize is None) or (self.getLoadLocation is None):
-                return
-        self.shaft.updateLoad(self.getFocusedColumn(), loadType = text)
 
-    def setLoadType(self, column, t):
+    def setConstraintType(self, column, t):
-        self.setListValue("LoadType", column, t)
+        self.setListValue("ConstraintType", column, t)
-        self.shaft.updateLoad(column, loadType = t)
+        self.shaft.updateConstraint(column, t)
 
-    def getLoadType(self, column):
+    def getConstraintType(self, column):
-        return self.getListValue("LoadType", column)
+        return self.getListValue("ConstraintType", column)
-
-    def setLoadSize(self, column, s):
-        self.setDoubleValue("LoadSize", column, s)
-        self.shaft.updateLoad(column, loadSize = s)
-
-    def getLoadSize(self, column):
-        return self.getDoubleValue("LoadSize", column)
-
-    def setLoadLocation(self, column, l):
-        self.setDoubleValue("LoadLocation", column, l)
-        self.shaft.updateLoad(column, loadLocation = l)
-
-    def getLoadLocation(self, column):
-        return self.getDoubleValue("LoadLocation", column)
 
     def slotStartEdgeType(self, old, new):
         pass
@@ -334,7 +309,7 @@ class WizardShaftTable:
     def getListValue(self, row, column):
         widget = self.widget.cellWidget(self.rowDict[row], column)
         if widget is not None:
-            return widget.currentText().toAscii()[0].upper()
+            return widget.currentText().toAscii() #[0].upper()
         else:
             return None