system: auto relax some redundant composit constraints

The order of appearance of the constraints determines how redundancies
are handled.

PlaneCoincident: the second PlaneCoincident of the same pair or parts is
relax to a 2d point horizontal or vertical constraint. The third and
onwards are ignored.  Any PlaneCoincident constraint will be ignored if
there are any existing PlaneCoincident or AxialAlignment.

PlaneAlignment: the second one is relaxed to two PointInPlane or
PointPlaneDistance constraint. The third one is relaxed to one
PointInPlane or PointPlaneDistance. The forth and onwards are ignored.
Any PlaneCoincident constraint will be ignored if there are any existing
PlaneCoincident or AxialAlignment.

AxialAlignment: ignored if there are any existing PlaneAlignment or
PlaneCoincident constraints

Author: realthunder

constraint.py

@@ -47,8 +47,9 @@ def _p(solver,partInfo,subname,shape,retAll=False):
         system.log('cache {}: {}'.format(key,h))
         return h if retAll else h[0]
 
+    v = utils.getElementPos(shape)
+
     if utils.isDraftWire(part):
-        v = utils.getElementPos(shape)
         nameTag = partInfo.PartName + '.' + key
         v = partInfo.Placement.multVec(v)
         params = []
@@ -83,7 +84,6 @@ def _p(solver,partInfo,subname,shape,retAll=False):
         system.log('{}: add circle point {},{}'.format(key,h,e))
 
     else:
-        v = utils.getElementPos(shape)
         nameTag = partInfo.PartName + '.' + key
         system.NameTag = nameTag
         e = system.addPoint3dV(*v)
@@ -92,6 +92,8 @@ def _p(solver,partInfo,subname,shape,retAll=False):
         h = [h,e]
         system.log('{}: {},{}'.format(key,h,partInfo.Group))
 
+    # use the point entity as key to store its original position vector
+    partInfo.EntityMap[h[0]] = [v]
     partInfo.EntityMap[key] = h
     return h if retAll else h[0]
 
@@ -130,7 +132,14 @@ def _n(solver,partInfo,subname,shape,retAll=False):
         system.NameTag = nameTag + 'xt'
         h.append(system.addTransform(e,*partInfo.Params,group=partInfo.Group))
 
+        # also add local x pointing vector (1,0,0) 
+        px = rot.multVec(FreeCAD.Vector(1,0,0))
+        hx = system.addPoint3dV(px.x,px.y,px.z)
+        h.append(system.addTransform(hx,*partInfo.Params,group=partInfo.Group))
+
         system.log('{}: {},{}'.format(key,h,partInfo.Group))
+        # use the normal entity as key to store its original rotation
+        partInfo.EntityMap[h[0]] = [rot]
         partInfo.EntityMap[key] = h
     return h if retAll else h[0]
 
@@ -235,9 +244,8 @@ def _w(solver,partInfo,subname,shape,retAll=False):
         n = _n(solver,partInfo,subname,shape,True)
         system.NameTag = partInfo.PartName + '.' + key
         w = system.addWorkplane(p,n[0],group=partInfo.Group)
-        h = [w,p] + n
+        h = [w,p,n]
         system.log('{}: {},{}'.format(key,h,partInfo.Group))
-        partInfo.EntityMap[key] = h
     return h if retAll else h[0]
 
 def _wa(solver,partInfo,subname,shape,retAll=False):
@@ -269,6 +277,7 @@ def _c(solver,partInfo,subname,shape,requireArc=False,retAll=False):
     if utils.isDraftCircle(partInfo.Part):
         part = partInfo.Part
         w,p,n = partInfo.Workplane[:3]
+        n = n[0]
 
         if system.sketchPlane and not solver.isFixedElement(part,subname):
             system.NameTag = nameTag + '.o'
@@ -313,6 +322,7 @@ def _c(solver,partInfo,subname,shape,requireArc=False,retAll=False):
             partInfo.EntityMap[sub] = h
     else:
         w,p,n = _w(solver,partInfo,subname,shape,True)[:3]
+        n = n[0]
         r = utils.getElementCircular(shape)
         if not r:
             raise RuntimeError('shape is not cicular')
@@ -831,14 +841,17 @@ def prepare(cls,obj,solver):
             return
         e0 = None
         ret = []
+        firstInfo = None
         for e in elements:
             info = e.Proxy.getInfo()
             partInfo = solver.getPartInfo(info)
             if not e0:
                 e0 = cls._entityDef[0](solver,partInfo,info.Subname,info.Shape)
+                firstInfo = partInfo
             else:
                 e = cls._entityDef[0](solver,partInfo,info.Subname,info.Shape)
                 params = props + [e0,e]
+                solver.system.checkRedundancy(obj,firstInfo,partInfo)
                 h = func(*params,group=solver.group)
                 if isinstance(h,(list,tuple)):
                     ret += list(h)
@@ -873,6 +886,7 @@ def prepare(cls,obj,solver):
                 params = props + [e1,e2]
             else:
                 params = props + [e2,e1]
+            solver.system.checkRedendancy(obj,prevInfo,partInfo)
             h = func(*params,group=solver.group)
             if isinstance(h,(list,tuple)):
                 ret += list(h)
@@ -892,7 +906,6 @@ class PlaneCoincident(BaseCascade):
         'Add a "{}" constraint to conincide planes of two or more parts.\n'\
         'The planes are coincided at their centers with an optional distance.'
 
-
 class PlaneAlignment(BaseCascade):
     _id = 37
     _iconName = 'Assembly_ConstraintAlignment.svg'

solver.py

@@ -17,8 +17,11 @@
 #            plane of the part.
 # EntityMap: string -> entity handle map, for caching
 # Group: transforming entity group handle
+# CstrMap: map from other part to the constrains between this and the othe part.
+#          This is for auto constraint DOF reduction. Only some composite
+#          constraints will be mapped.
 PartInfo = namedtuple('SolverPartInfo', ('Part','PartName','Placement',
-    'Params','Workplane','EntityMap','Group'))
+    'Params','Workplane','EntityMap','Group','CstrMap'))
 
 class Solver(object):
     def __init__(self,assembly,reportFailed,dragPart,recompute,rollback):
@@ -37,13 +40,17 @@ def __init__(self,assembly,reportFailed,dragPart,recompute,rollback):
 
         self.system.GroupHandle = self._fixedGroup
 
-        # convenience constant of zero
+        # convenience constant of zero and one
         self.v0 = self.system.addParamV(0,group=self._fixedGroup)
+        self.v1 = self.system.addParamV(1,group=self._fixedGroup)
 
         # convenience normals
         rotx = FreeCAD.Rotation(FreeCAD.Vector(0,1,0),-90)
         self.nx = self.system.addNormal3dV(*utils.getNormal(rotx))
 
+        # convenience x pointing vector
+        self.px = self.system.addPoint3d(self.v1,self.v0,self.v0)
+
         self._fixedParts = Constraint.getFixedParts(self,cstrs)
         for part in self._fixedParts:
             self._fixedElements.add((part,None))
@@ -55,7 +62,8 @@ def __init__(self,assembly,reportFailed,dragPart,recompute,rollback):
             if ret:
                 if isinstance(ret,(list,tuple)):
                     for h in ret:
-                        self._cstrMap[h] = cstr
+                        if not isinstance(h,(list,tuple)):
+                            self._cstrMap[h] = cstr
                 else:
                     self._cstrMap[ret] = cstr
 
@@ -217,17 +225,20 @@ def getPartInfo(self,info,fixed=False,group=0):
             self.system.NameTag = info.PartName + '.nx'
             nx = self.system.addTransform(self.nx,
                     self.v0,self.v0,self.v0,*params[3:],group=g)
+            px = self.system.addTransform(self.px,self.v0,self.v0,
+                    self.v0,*params[3:],group=g)
             self.system.NameTag = info.PartName + '.w'
             w = self.system.addWorkplane(p,n,group=g)
-            h = (w,p,n,nx)
+            h = (w,p,(n,nx,px))
 
         partInfo = PartInfo(Part = info.Part,
                             PartName = info.PartName,
                             Placement = info.Placement.copy(),
                             Params = params,
                             Workplane = h,
                             EntityMap = {},
-                            Group = group if group else g)
+                            Group = group if group else g,
+                            CstrMap = {})
 
         self.system.log('{}, {}'.format(partInfo,g))
 

system.py

@@ -1,5 +1,6 @@
 import os
-from .utils import getIcon, syslogger as logger, objName
+from .constraint import cstrName
+from .utils import getIcon, syslogger as logger, objName, project2D
 from .proxy import ProxyType, PropertyInfo
 
 class System(ProxyType):
@@ -58,12 +59,12 @@ def getSystem(mcs,obj):
             return proxy.getSystem(obj)
 
     @classmethod
-    def isConstraintSupported(mcs,obj,cstrName):
-        if cstrName == 'Locked':
+    def isConstraintSupported(mcs,obj,name):
+        if name == 'Locked':
             return True
         proxy = mcs.getProxy(obj)
         if proxy:
-            return proxy.isConstraintSupported(cstrName)
+            return proxy.isConstraintSupported(name)
 
 def _makePropInfo(name,tp,doc=''):
     PropertyInfo(System,name,tp,doc,group='Solver')
@@ -112,6 +113,12 @@ def __init__(self):
         super(SystemExtension,self).__init__()
         self.NameTag = ''
         self.sketchPlane = None
+        self.cstrObj = None
+        self.firstInfo = None
+        self.secondInfo = None
+
+    def checkRedundancy(self,obj,firstInfo,secondInfo):
+        self.cstrObj,self.firstInfo,self.secondInfo=obj,firstInfo,secondInfo
 
     def addSketchPlane(self,*args,**kargs):
         _ = kargs
@@ -127,12 +134,71 @@ def setOrientation(self,h,lockAngle,angle,n1,n2,nx1,nx2,group):
                 h.append(self.addAngle(angle,False,nx1,nx2,group=group))
         return h
 
+    def reportRedundancy(self,warn=False):
+        msg = '{} between {} and {}'.format(cstrName(self.cstrObj),
+                self.firstInfo.PartName, self.secondInfo.PartName)
+        if warn:
+            logger.warn('skip redundant ' + msg)
+        else:
+            logger.info('auto relax ' + msg)
+
+    def countConstraints(self,increment,count,*names):
+        first,second = self.firstInfo,self.secondInfo
+        if not first or not second:
+            return
+        ret = 0
+        for name in names:
+            cstrs = first.CstrMap.get(second.Part,{}).get(name,None)
+            if not cstrs:
+                if increment:
+                    cstrs = second.CstrMap.setdefault(
+                                first.Part,{}).setdefault(name,[])
+                else:
+                    cstrs = second.CstrMap.get(first.Part,{}).get(name,[])
+            if increment:
+                cstrs += [None]*increment
+            ret += len(cstrs)
+            if ret >= count:
+                if ret>count:
+                    self.reportRedundancy(True)
+                    return -1
+                else:
+                    self.reportRedundancy()
+        return ret
+
     def addPlaneCoincident(self,d,lockAngle,angle,e1,e2,group=0):
         if not group:
             group = self.GroupHandle
-        w1,p1,n1,nx1 = e1[:4]
-        _,p2,n2,nx2 = e2[:4]
+        w1,p1,n1 = e1[:3]
+        _,p2,n2 = e2[:3]
+        n1,nx1 = n1[:2]
+        n2,nx2 = n2[:2]
         h = []
+        count = self.countConstraints(0,1,'Alignment','Axial')
+        if count<0:
+            # We do not allow plane coincident coexist with other composite
+            # constraints
+            return
+        count = self.countConstraints(2 if lockAngle else 1,2,'Coincident')
+        if count<0:
+            return
+        if not lockAngle and count==2:
+            # if there is already some other plane coincident constraint set for
+            # this pair of parts, we reduce this second constraint to either a
+            # points horizontal or vertical constraint, i.e. reduce the
+            # constraining DOF down to 1.
+            #
+            # We project the initial points to the first element plane, and
+            # check for differences in x and y components of the points to
+            # determine whether to use horizontal or vertical constraint.
+            v1,v2 = project2D(self.firstInfo.EntityMap[n1][0],
+                              self.firstInfo.EntityMap[p1][0],
+                              self.secondInfo.EntityMap[p2][0])
+            if abs(v1.x-v2.x) < abs(v1.y-v2.y):
+                h.append(self.addPointsHorizontal(p1,p2,w1,group=group))
+            else:
+                h.append(self.addPointsVertical(p1,p2,w1,group=group))
+            return h
         if d:
             h.append(self.addPointPlaneDistance(d,p2,w1,group=group))
             h.append(self.addPointsCoincident(p1,p2,w1,group=group))
@@ -143,20 +209,41 @@ def addPlaneCoincident(self,d,lockAngle,angle,e1,e2,group=0):
     def addPlaneAlignment(self,d,lockAngle,angle,e1,e2,group=0):
         if not group:
             group = self.GroupHandle
-        w1,_,n1,nx1 = e1[:4]
-        _,p2,n2,nx2 = e2[:4]
+        w1,_,n1 = e1[:4]
+        _,p2,n2 = e2[:4]
+        n1,nx1 = n1[:2]
+        n2,nx2,px2 = n2[:3]
         h = []
+        count = self.countConstraints(0,1,'Coincident','Axial')
+        if count<0:
+            return
+        count = self.countConstraints(2 if lockAngle else 1,3,'Alignment')
+        if count<0:
+            return
         if d:
             h.append(self.addPointPlaneDistance(d,p2,w1,group=group))
         else:
             h.append(self.addPointInPlane(p2,w1,group=group))
-        return self.setOrientation(h,lockAngle,angle,n1,n2,nx1,nx2,group)
+        if count==1 or (lockAngle and count==2):
+            h.append(self.setOrientation(h,lockAngle,angle,n1,n2,nx1,nx2,group))
+        elif count==2:
+            self.reportRedundancy()
+            if d:
+                h.append(self.addPointPlaneDistance(d,px2,w1,group=group))
+            else:
+                h.append(self.addPointInPlane(px2,w1,group=group))
+        return h
 
     def addAxialAlignment(self,lockAngle,angle,e1,e2,group=0):
         if not group:
             group = self.GroupHandle
-        w1,p1,n1,nx1 = e1[:4]
-        _,p2,n2,nx2 = e2[:4]
+        count = self.countConstraints(0,1,'Coincident','Alignment')
+        if count<0:
+            return
+        w1,p1,n1 = e1[:3]
+        _,p2,n2 = e2[:3]
+        n1,nx1 = n1[:2]
+        n2,nx2 = n2[:2]
         h = []
         h.append(self.addPointsCoincident(p1,p2,w1,group=group))
         return self.setOrientation(h,lockAngle,angle,n1,n2,nx1,nx2,group)

utils.py

@@ -572,3 +572,9 @@ def getLabel(obj):
                 label = label[:-i]
             break
     return label
+
+def project2D(rot,*vectors):
+    vx = rot.multVec(FreeCAD.Vector(1,0,0))
+    vy = rot.multVec(FreeCAD.Vector(0,1,0))
+    return [FreeCAD.Vector(v.dot(vx),v.dot(vy),0) for v in vectors]
+