diff --git a/FCADLogger.py b/FCADLogger.py
new file mode 100644
index 0000000..e2edbcb
--- /dev/null
+++ b/FCADLogger.py
@@ -0,0 +1,66 @@
+import os
+from datetime import datetime
+import inspect
+import FreeCAD, FreeCADGui
+
+class FCADLogger:
+    def __init__(self, tag, **kargs):
+        self.tag = tag
+        self.levels = { 'error':0, 'warn':1, 'info':2,
+                'debug':3, 'trace':4 }
+        self.printer = [
+                FreeCAD.Console.PrintError,
+                FreeCAD.Console.PrintWarning,
+                FreeCAD.Console.PrintMessage,
+                FreeCAD.Console.PrintLog,
+                FreeCAD.Console.PrintLog ]
+        self.laststamp = datetime.now()
+        for key in ('printTag','updateUI','timing','lineno'):
+            setattr(self,key,kargs.get(key,True))
+
+    def _isEnabledFor(self,level):
+        return FreeCAD.getLogLevel(self.tag) >= level
+
+    def isEnabledFor(self,level):
+        self._isEnabledOf(self.levels[level])
+
+    def error(self,msg,frame=0):
+        self.log(0,msg,frame+1)
+
+    def warn(self,msg,frame=0):
+        self.log(1,msg,frame+1)
+
+    def info(self,msg,frame=0):
+        self.log(2,msg,frame+1)
+
+    def debug(self,msg,frame=0):
+        self.log(3,msg,frame+1)
+
+    def trace(self,msg,frame=0):
+        self.log(4,msg,frame+1)
+
+    def log(self,level,msg,frame=0):
+        if not self._isEnabledFor(level):
+            return
+
+        prefix = ''
+
+        if self.printTag:
+            prefix += '<{}> '.format(self.tag)
+
+        if self.timing:
+            now = datetime.now()
+            prefix += '{} - '.format((now-self.laststamp).total_seconds())
+            self.laststamp = now
+
+        if self.lineno:
+            stack = inspect.stack()[frame+1]
+            prefix += '{}({}): '.format(os.path.basename(stack[1]),stack[2])
+
+        self.printer[level]('{}{}\n'.format(prefix,msg))
+
+        if self.updateUI:
+            try:
+                FreeCADGui.updateGui()
+            except Exception:
+                pass
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..61c20a0
--- /dev/null
+++ b/__init__.py
@@ -0,0 +1,26 @@
+
+import FreeCAD, FreeCADGui, Part
+import asm3.assembly as assembly
+import asm3.constraint as constraint
+import asm3.utils as utils
+import asm3.solver as solver
+from asm3.assembly import Assembly,AsmConstraint
+
+def test():
+    doc = FreeCAD.newDocument()
+    cylinder1 = doc.addObject('Part::Cylinder','cylinder1')
+    cylinder1.Visibility = False
+    asm1 = Assembly.make(doc)
+    asm1.Proxy.getPartGroup().setLink({-1:cylinder1})
+    cylinder2 = doc.addObject('Part::Cylinder','cylinder2')
+    cylinder2.Visibility = False
+    asm2 = Assembly.make(doc)
+    asm2.Placement.Base.z = -20
+    asm2.Proxy.getPartGroup().setLink({-1:cylinder2})
+    doc.recompute()
+    FreeCADGui.SendMsgToActiveView("ViewFit")
+    asm = Assembly.make(doc)
+    asm.Proxy.getPartGroup().setLink((asm1,asm2))
+    asm1.Visibility = False
+    asm2.Visibility = False
+
diff --git a/assembly.py b/assembly.py
new file mode 100644
index 0000000..c664e8b
--- /dev/null
+++ b/assembly.py
@@ -0,0 +1,982 @@
+import sys, os
+sys.path.append(os.path.dirname(os.path.realpath(__file__)))
+
+from collections import namedtuple
+
+import FreeCAD, FreeCADGui
+
+import asm3.constraint as constraint
+from asm3.utils import logger, objName
+
+def setupUndo(doc,undoDocs,name='Assembly3 solve'):
+    if doc in undoDocs:
+        return
+    doc.openTransaction(name)
+    undoDocs.add(doc)
+
+def isTypeOf(obj,tp,resolve=False):
+    if not obj:
+        return False
+    if not tp:
+        return True
+    if resolve:
+        obj = obj.getLinkedObject(True)
+    return isinstance(getattr(obj,'Proxy',None),tp)
+
+def checkType(obj,tp,resolve=False):
+    if not isTypeOf(obj,tp,resolve):
+        raise TypeError('Expect object "{}" to be of type "{}"'.format(
+                objName(obj),tp.__name__))
+
+def getProxy(obj,tp):
+    checkType(obj,tp)
+    return obj.Proxy
+
+class AsmBase(object):
+    def __init__(self):
+        self.obj = None
+
+    def __getstate__(self):
+        return
+
+    def __setstate__(self,_state):
+        return
+
+    def attach(self,obj):
+        obj.addExtension('App::LinkBaseExtensionPython', None)
+        self.linkSetup(obj)
+
+    def linkSetup(self,obj):
+        assert getattr(obj,'Proxy',None)==self
+        self.obj = obj
+        return
+
+    def getViewProviderName(self,_obj):
+        return 'Gui::ViewProviderLinkPython'
+
+    def onDocumentRestored(self, obj):
+        self.linkSetup(obj)
+
+    def onChanged(self,_obj,_prop):
+        pass
+
+class ViewProviderAsmBase(object):
+    def __init__(self,vobj):
+        vobj.Visibility = False
+        vobj.Proxy = self
+        self.attach(vobj)
+
+    def attach(self,vobj):
+        self.ViewObject = vobj
+
+    def __getstate__(self):
+        return None
+
+    def __setstate__(self, _state):
+        return None
+
+
+class AsmGroup(AsmBase):
+    def linkSetup(self,obj):
+        super(AsmGroup,self).linkSetup(obj)
+        obj.configLinkProperty(
+                'VisibilityList',LinkMode='GroupMode',ElementList='Group')
+        self.setGroupMode()
+
+    def setGroupMode(self):
+        self.obj.GroupMode = 1 # auto delete children
+        self.obj.setPropertyStatus('GroupMode','Hidden')
+        self.obj.setPropertyStatus('GroupMode','Immutable')
+        self.obj.setPropertyStatus('GroupMode','Transient')
+
+    def attach(self,obj):
+        obj.addProperty("App::PropertyLinkList","Group","Base",'')
+        obj.addProperty("App::PropertyBoolList","VisibilityList","Base",'')
+        obj.addProperty("App::PropertyEnumeration","GroupMode","Base",'')
+        super(AsmGroup,self).attach(obj)
+
+
+class ViewProviderAsmGroup(ViewProviderAsmBase):
+
+    def claimChildren(self):
+        return self.ViewObject.Object.Group
+
+
+class AsmPartGroup(AsmGroup):
+    def __init__(self,parent):
+        self.parent = getProxy(parent,Assembly)
+        super(AsmPartGroup,self).__init__()
+
+    def setGroupMode(self):
+        pass
+
+    @staticmethod
+    def make(parent,name='Parts'):
+        obj = parent.Document.addObject("App::FeaturePython",name,
+                    AsmPartGroup(parent),None,True)
+        ViewProviderAsmPartGroup(obj.ViewObject)
+        return obj
+
+
+class ViewProviderAsmPartGroup(ViewProviderAsmBase):
+    def onDelete(self,_obj,_subs):
+        return False
+
+
+class AsmElement(AsmBase):
+    def __init__(self,parent):
+        self.shape = None
+        self.parent = getProxy(parent,AsmElementGroup)
+        super(AsmElement,self).__init__()
+
+    def linkSetup(self,obj):
+        super(AsmElement,self).linkSetup(obj)
+        obj.configLinkProperty('LinkedObject')
+        obj.setPropertyStatus('LinkedObject','Immutable')
+        obj.setPropertyStatus('LinkedObject','ReadOnly')
+
+    def attach(self,obj):
+        obj.addProperty("App::PropertyXLink","LinkedObject"," Link",'')
+        super(AsmElement,self).attach(obj)
+
+    def execute(self,_obj):
+        self.getShape(True)
+        return False
+
+    def getShape(self,refresh=False):
+        if not refresh:
+            ret = getattr(self,'shape',None)
+            if ret:
+                return ret
+        self.shape = None
+        self.shape = self.obj.getSubObject('')
+        return self.shape
+
+    def getAssembly(self):
+        return self.parent.parent
+
+    def getSubElement(self):
+        link = self.obj.LinkedObject
+        if isinstance(link,tuple):
+            return link[1].split('.')[-1]
+        return ''
+
+    def getSubName(self):
+        link = self.obj.LinkedObject
+        if not isinstance(link,tuple):
+            raise RuntimeError('Invalid element link "{}"'.format(
+                objName(self.obj)))
+        return link[1]
+
+    def setLink(self,owner,subname):
+        # subname must be relative to the part group object of the parent
+        # assembly
+
+        # check old linked object for auto re-label
+        obj = self.obj
+        linked = obj.getLinkedObject(False)
+        if linked and linked!=obj:
+            label = linked.Label + '_' + self.getSubElement()
+        else:
+            label = ''
+
+        obj.setLink(owner,subname)
+
+        if obj.Label==obj.Name or obj.Label==label:
+            linked = obj.getLinkedObject(False)
+            if linked and linked!=obj:
+                obj.Label = linked.Label+'_'+self.getSubElement()
+            else:
+                obj.Label = obj.Name
+
+    Selection = namedtuple('AsmElementSelection',
+            ('Assembly','Element','Subname'))
+
+    @staticmethod
+    def getSelection():
+        '''
+        Parse Gui.Selection for making a element
+
+        If there is only one selection, then the selection must refer to a sub
+        element of some part object of an assembly. We shall create a new
+        element beloning to the top-level assembly
+
+        If there are two selections, then first one shall be either the
+        element group or an individual element. The second selection shall
+        be a sub element belong to a child assembly of the parent assembly of
+        the first selected element/element group
+        '''
+        sels = FreeCADGui.Selection.getSelectionEx('',False)
+        if not sels:
+            return
+        if len(sels)>1:
+            raise RuntimeError(
+                    'The selections must have a common (grand)parent assembly')
+
+        sel = sels[0]
+        subs = sel.SubElementNames
+        if len(subs)>2:
+            raise RuntimeError('At most two selection is allowed.\n'
+                'The first selection must be a sub element belonging to some '
+                'assembly. The optional second selection must be an element '
+                'belonging to the same assembly of the first selection')
+
+        subElement = subs[0].split('.')[-1]
+        if not subElement:
+            raise RuntimeError(
+                'Please select a sub element belonging to some assembly')
+
+        link = Assembly.findPartGroup(sel.Object,subs[0])
+        if not link:
+            raise RuntimeError(
+                    'Selected sub element does not belong to an assembly')
+
+        element = None
+        if len(subs)>1:
+            ret = Assembly.findElementGroup(sel.Object,subs[1])
+            if not ret:
+                raise RuntimeError('The second selection must be an element')
+
+            if ret.Assembly != link.Assembly:
+                raise RuntimeError(
+                        'The two selections must belong to the same assembly')
+
+            element = ret.Object.getSubObject(ret.Subname,1)
+            if not isTypeOf(element,AsmElement):
+                raise RuntimeError('The second selection must be an element')
+
+        return AsmElement.Selection(
+                link.Assembly,element,link.Subname+subElement)
+
+    @staticmethod
+    def make(selection=None,name='Element'):
+        if not selection:
+            selection = AsmElement.getSelection()
+        assembly = getProxy(selection.Assembly,Assembly)
+        element = selection.Element
+        if not element:
+            elements = assembly.getElementGroup()
+            # try to search the element group for an existing element
+            for e in elements.Group:
+                if getProxy(e,AsmElement).getSubName() == selection.Subname:
+                    return element
+            element = elements.Document.addObject("App::FeaturePython",
+                    name,AsmElement(elements),None,True)
+            ViewProviderAsmElement(element.ViewObject)
+            elements.setLink({-1:element})
+
+        getProxy(element,AsmElement).setLink(
+                assembly.getPartGroup(),selection.Subname)
+        return element
+
+
+class ViewProviderAsmElement(ViewProviderAsmBase):
+    def attach(self,vobj):
+        super(ViewProviderAsmElement,self).attach(vobj)
+        vobj.OverrideMaterial = True
+        vobj.ShapeMaterial.DiffuseColor = self.getDefaultColor()
+        vobj.ShapeMaterial.EmissiveColor = self.getDefaultColor()
+        vobj.DrawStyle = 1
+        vobj.LineWidth = 4
+        vobj.PointSize = 6
+
+    def getDefaultColor(self):
+        return (60.0/255.0,1.0,1.0)
+
+
+class AsmElementLink(AsmBase):
+    def __init__(self,parent):
+        super(AsmElementLink,self).__init__()
+        self.info = None
+        self.parent = getProxy(parent,AsmConstraint)
+
+    def linkSetup(self,obj):
+        super(AsmElementLink,self).linkSetup(obj)
+        obj.configLinkProperty('LinkedObject')
+        obj.setPropertyStatus('LinkedObject','Immutable')
+        obj.setPropertyStatus('LinkedObject','ReadOnly')
+
+    def attach(self,obj):
+        obj.addProperty("App::PropertyXLink","LinkedObject"," Link",'')
+        super(AsmElementLink,self).attach(obj)
+
+    def execute(self,_obj):
+        self.getInfo(True)
+        return False
+
+    def getAssembly(self):
+        return self.parent.parent.parent
+
+    def getElement(self):
+        linked = self.obj.getLinkedObject(False)
+        if not linked:
+            raise RuntimeError('Element link broken')
+        if not isTypeOf(linked,AsmElement):
+            raise RuntimeError('Invalid element type')
+        return linked.Proxy
+
+    def getSubName(self):
+        link = self.obj.LinkedObject
+        if not isinstance(link,tuple):
+            raise RuntimeError('Invalid element link "{}"'.format(
+                objName(self.obj)))
+        return link[1]
+
+    def getShapeSubName(self):
+        element = self.getElement()
+        assembly = element.getAssembly()
+        if assembly == self.getAssembly():
+            return element.getSubName()
+        # pop two names from back (i.e. element group, element)
+        subname = self.getSubName()
+        sub = subname.split('.')[:-3]
+        sub = '.'.join(sub) + '.' + assembly.getPartGroup().Name + \
+              '.' + element.getSubName()
+        logger.debug('shape subname {} -> {}'.format(subname,sub))
+        return sub
+
+    def prepareLink(self,owner,subname):
+        assembly = self.getAssembly()
+        sobj = owner.getSubObject(subname,1)
+        if not sobj:
+            raise RuntimeError('invalid element link {} broken: {}'.format(
+                objName(owner),subname))
+        if isTypeOf(sobj,AsmElementLink):
+            # if it points to another AsElementLink that belongs the same
+            # assembly, simply return the same link
+            if sobj.Proxy.getAssembly() == assembly:
+                return (owner,subname)
+            # If it is from another assembly (i.e. a nested assembly), convert
+            # the subname reference by poping three names (constraint group,
+            # constraint, element link) from the back, and then append with the
+            # element link's own subname reference
+            sub = subname.split('.')[:-4]
+            sub = '.'.join(subname)+'.'+sobj.Proxy.getSubName()
+            logger.debug('convert element link {} -> {}'.format(subname,sub))
+            return (owner,sub)
+
+        if isTypeOf(sobj,AsmElement):
+            return (owner,subname)
+
+        # try to see if the reference comes from some nested assembly
+        ret = assembly.findChild(owner,subname,recursive=True)
+        if not ret:
+            # It is from a non assembly child part, then use our own element
+            # group as the holder for elements
+            ret = [Assembly.Selection(assembly.obj,owner,subname)]
+
+        if not isTypeOf(ret[-1].Object,AsmPartGroup):
+            raise RuntimeError('Invalid element link ' + subname)
+
+        # call AsmElement.make to either create a new element, or an existing
+        # element if there is one
+        element = AsmElement.make(AsmElement.Selection(
+                    ret[-1].Assembly,None,ret[-1].Subname))
+        if ret[-1].Assembly == assembly.obj:
+            return (assembly.getElementGroup(),element.Name+'.')
+
+        elementSub = ret[-1].Object.Name + '.' + ret[-1].Subname
+        sub = subname[:-(len(elementSub)+1)] + '.' + \
+            ret[-1].Assembly.Proxy.getElementGroup().Name + '.' + \
+            element.Name + '.'
+        logger.debug('generate new element {} -> {}'.format(subname,sub))
+        return (owner,sub)
+
+    def setLink(self,owner,subname):
+        obj = self.obj
+        obj.setLink(*self.prepareLink(owner,subname))
+        linked = obj.getLinkedObject(False)
+        if linked and linked!=obj:
+            obj.Label = 'Link_'+linked.Label
+        else:
+            obj.Label = obj.Name
+
+    Info = namedtuple('AsmElementLinkInfo',
+            ('Part','PartName','Placement','Object','Subname','Shape'))
+
+    def getInfo(self,refresh=False):
+        if not refresh:
+            ret = getattr(self,'info',None)
+            if ret:
+                return ret
+        self.info = None
+        assembly = self.getAssembly()
+        subname = self.getShapeSubName()
+        names = subname.split('.')
+        partGroup = assembly.getPartGroup()
+
+        part = partGroup.getSubObject(names[0]+'.',1)
+        if not part:
+            raise RuntimeError('Eelement link "{}" borken: {}'.format(
+                objName(self.obj),subname))
+
+        # For storing the shape of the element with proper transformation
+        shape = None
+        # For storing the placement of the movable part
+        pla = None
+        # For storing the actual geometry object of the part, in case 'part' is
+        # a link
+        obj = None
+
+        if not isTypeOf(part,Assembly,True) and part!=partGroup.Group[0]:
+            getter = getattr(part.getLinkedObject(True),'getLinkExtProperty')
+
+            # special treatment of link array (i.e. when ElementCount!=0), we
+            # allow the array element to be moveable by the solver
+            if getter and getter('ElementCount'):
+
+                # store both the part (i.e. the link array), and the array
+                # element object
+                part = (part,part.getSubObject(names[1]+'.',1))
+
+                # trim the subname to be after the array element
+                sub = '.'.join(names[2:])
+
+                # There are two states of an link array. 
+                if getter('ElementList'):
+                    # a) The elements are expanded as individual objects, i.e
+                    # when ElementList has members, then the moveable Placement
+                    # is a property of the array element. So we obtain the shape
+                    # before 'Placement' by setting 'transform' set to False.
+                    shape=part[1].getSubObject(sub,transform=False)
+                    pla = part[1].Placement
+                    obj = part[0].getLinkedObject(False)
+                    partName = part[1].Name
+                else:
+                    # b) The elements are collapsed. Then the moveable Placement
+                    # is stored inside link object's PlacementList property. So,
+                    # the shape obtained below is already before 'Placement',
+                    # i.e. no need to set 'transform' to False.
+                    shape=part[1].getSubObject(sub)
+                    obj = part[1]
+                    try:
+                        idx = names[1].split('_i')[-1]
+                        # we store the array index instead, in order to modified
+                        # Placement later when the solver is done. Also because
+                        # that when the elements are collapsed, there is really
+                        # no element object here.
+                        part = (part[0],int(idx),part[1])
+                        pla = part[0].PlacementList[idx]
+                    except ValueError:
+                        raise RuntimeError('invalid array subname of element '
+                            '"{}": {}'.format(objName(self.obj),subname))
+
+                    partName = '{}.{}.'.format(part[0].Name,idx)
+
+                subname = sub
+
+        if not shape:
+            # Here means, either the 'part' is an assembly or it is a non array
+            # object. We trim the subname reference to be after the part object.
+            # And obtain the shape before part's Placement by setting
+            # 'transform' to False
+            subname = '.'.join(names[1:])
+            shape = part.getSubObject(subname,transform=False)
+            pla = part.Placement
+            obj = part.getLinkedObject(False)
+            partName = part.Name
+
+        self.info = AsmElementLink.Info(
+                part,partName,pla.copy(),obj,subname,shape.copy())
+        return self.info
+
+    @staticmethod
+    def setPlacement(part,pla,undoDocs):
+        '''
+        called by solver after solving to adjust the placement.
+        
+        part: obtained by AsmConstraint.getInfo().Part
+        pla: the new placement
+        '''
+        if isinstance(part,tuple):
+            if isinstance(part[1],int):
+                setupUndo(part[0].Document,undoDocs)
+                part[0].PlacementList = {part[1]:pla}
+            else:
+                setupUndo(part[1].Document,undoDocs)
+                part[1].Placement = pla
+        else:
+            setupUndo(part.Document,undoDocs)
+            part.Placement = pla
+
+    MakeInfo = namedtuple('AsmElementLinkSelection',
+            ('Constraint','Owner','Subname'))
+
+    @staticmethod
+    def make(info,name='ElementLink'):
+        element = info.Constraint.Document.addObject("App::FeaturePython",
+                    name,AsmElementLink(info.Constraint),None,True)
+        ViewProviderAsmElementLink(element.ViewObject)
+        info.Constraint.setLink({-1:element})
+        element.Proxy.setLink(info.Owner,info.Subname)
+        return element
+
+
+class ViewProviderAsmElementLink(ViewProviderAsmBase):
+    pass
+
+
+class AsmConstraint(AsmGroup):
+
+    def __init__(self,parent):
+        self.elements = None
+        self.parent = getProxy(parent,AsmConstraintGroup)
+        super(AsmConstraint,self).__init__()
+
+    def attach(self,obj):
+        # Property '_Type' is hidden from editor. The type is for the solver to 
+        # store some internal type id of the constraint, to avoid potential
+        # problem of version upgrade in the future. The type id is oqaque to the
+        # objects in this module. The solve is reponsible to add the actual
+        # 'Type' enumeration property that is avaiable for user to change in the
+        # editor
+        obj.addProperty("App::PropertyInteger","_Type","Base",'',0,False,True)
+        super(AsmConstraint,self).attach(obj)
+
+    def onChanged(self,obj,prop):
+        constraint.onChanged(obj,prop)
+        super(AsmConstraint,self).onChanged(obj,prop)
+
+    def linkSetup(self,obj):
+        self.elements = None
+        super(AsmConstraint,self).linkSetup(obj)
+        obj.setPropertyStatus('VisibilityList','Output')
+        for o in obj.Group:
+            getProxy(o,AsmElementLink).parent = self
+        constraint.attach(obj)
+
+    def execute(self,_obj):
+        self.getElements(True)
+        return False
+
+    def getElements(self,refresh=False):
+        if refresh:
+            self.elements = None
+        ret = getattr(self,'elements',None)
+        obj = self.obj
+        if ret or not obj._Type:
+            return ret
+        shapes = []
+        elements = []
+        for o in obj.Group:
+            checkType(o,AsmElementLink)
+            info = o.Proxy.getInfo()
+            shapes.append(info.Shape)
+            elements.append(o)
+        constraint.check(obj._Type,shapes)
+        self.elements = elements
+        return self.elements
+
+    Selection = namedtuple('ConstraintSelection',
+                    ('Assembly','Constraint','Elements'))
+
+    @staticmethod
+    def getSelection(tp=0):
+        '''
+        Parse Gui.Selection for making a constraint
+
+        The selected elements must all belong to the same immediate parent
+        assembly. 
+        '''
+        sels = FreeCADGui.Selection.getSelectionEx('',False)
+        if not sels:
+            return
+        if len(sels)>1:
+            raise RuntimeError(
+                    'The selections must have a common (grand)parent assembly')
+
+        sel = sels[0]
+        cstr = None
+        elements = []
+        assembly = None
+        for sub in sel.SubElementNames:
+            sobj = sel.Object.getSubObject(sub,1)
+            ret = Assembly.findChild(sel.Object,sub,recursive=True)
+            if not ret:
+                raise RuntimeError('Selection {}.{} is not from an '
+                    'assembly'.format(sel.Object.Name,sub))
+            if not assembly:
+                # check if the selection is a constraint group or a constraint
+                if isTypeOf(sobj,AsmConstraintGroup):
+                    assembly = ret[-1].Assembly
+                    continue
+                if isTypeOf(sobj,AsmConstraint):
+                    cstr = sobj
+                    assembly = ret[-1].Assembly
+                    continue
+                assembly = ret[0].Assembly
+
+            found = None
+            for r in ret:
+                if r.Assembly == assembly:
+                    found = r
+                    break
+            if not found:
+                raise RuntimeError('Selection {}.{} is not from the target '
+                    'assembly {}'.format(sel.Object.Name,sub,objName(assembly)))
+
+            elements.append((found.Object,found.Subname))
+
+        check = None
+        if cstr and cstr._Type:
+            tp = cstr._Type
+            info = cstr.Proxy.getInfo()
+            check = [o.getShape() for o in info.Elements] + elements
+        elif tp:
+            check = elements
+        if check:
+            constraint.check(tp,check)
+
+        return AsmConstraint.Selection(assembly,cstr,elements)
+
+    @staticmethod
+    def make(tp=0, selection=None, name='Constraint'):
+        if not selection:
+            selection = AsmConstraint.getSelection(tp)
+        if selection.Constraint:
+            cstr = selection.Constraint
+        else:
+            constraints = selection.Assembly.Proxy.getConstraintGroup()
+            cstr = constraints.Document.addObject("App::FeaturePython",
+                    name,AsmConstraint(constraints),None,True)
+            ViewProviderAsmConstraint(cstr.ViewObject)
+            constraints.setLink({-1:cstr})
+            cstr._Type = tp
+
+        for e in selection.Elements:
+            AsmElementLink.make(AsmElementLink.MakeInfo(cstr,*e))
+        return cstr
+
+
+class ViewProviderAsmConstraint(ViewProviderAsmGroup):
+    def attach(self,vobj):
+        super(ViewProviderAsmConstraint,self).attach(vobj)
+        vobj.OverrideMaterial = True
+        vobj.ShapeMaterial.DiffuseColor = self.getDefaultColor()
+        vobj.ShapeMaterial.EmissiveColor = self.getDefaultColor()
+
+    def getDefaultColor(self):
+        return (1.0,60.0/255.0,60.0/255.0)
+
+
+class AsmConstraintGroup(AsmGroup):
+    def __init__(self,parent):
+        self.parent = getProxy(parent,Assembly)
+        super(AsmConstraintGroup,self).__init__()
+
+    def linkSetup(self,obj):
+        super(AsmConstraintGroup,self).linkSetup(obj)
+        obj.setPropertyStatus('VisibilityList','Output')
+        for o in obj.Group:
+            getProxy(o,AsmConstraint).parent = self
+
+    @staticmethod
+    def make(parent,name='Constraints'):
+        obj = parent.Document.addObject("App::FeaturePython",name,
+                AsmConstraintGroup(parent),None,True)
+        ViewProviderAsmConstraintGroup(obj.ViewObject)
+        return obj
+
+
+class ViewProviderAsmConstraintGroup(ViewProviderAsmBase):
+    pass
+
+
+class AsmElementGroup(AsmGroup):
+    def __init__(self,parent):
+        self.parent = getProxy(parent,Assembly)
+        super(AsmElementGroup,self).__init__()
+
+    def linkSetup(self,obj):
+        super(AsmElementGroup,self).linkSetup(obj)
+        obj.setPropertyStatus('VisibilityList','Output')
+        for o in obj.Group:
+            getProxy(o,AsmElement).parent = self
+
+    @staticmethod
+    def make(parent,name='Elements'):
+        obj = parent.Document.addObject("App::FeaturePython",name,
+                        AsmElementGroup(parent),None,True)
+        ViewProviderAsmElementGroup(obj.ViewObject)
+        return obj
+
+
+class ViewProviderAsmElementGroup(ViewProviderAsmBase):
+
+    def onDelete(self,_obj,_subs):
+        return False
+
+    def canDragObject(self,_obj):
+        return False
+
+    def canDragObjects(self):
+        return False
+
+    def canDragAndDropObject(self,_obj):
+        return False
+
+    def canDropObjectEx(self,_obj,owner,subname):
+        # check if is dropping a sub-element
+        if subname.rfind('.')+1 == len(subname):
+            return False
+        return self.ViewObject.Object.Proxy.parent.getPartGroup()==owner
+
+    def dropObjectEx(self,vobj,_obj,_owner,subname):
+        AsmElement.make(AsmElement.Selection(
+            vobj.Object.Proxy.parent.obj,None,subname))
+
+
+BuildShapeNames = ('No','Compound','Fuse','Cut')
+BuildShapeEnum = namedtuple('AsmBuildShapeEnum',BuildShapeNames)(
+        *range(len(BuildShapeNames)))
+
+
+class Assembly(AsmGroup):
+    def __init__(self):
+        self.constraints = None
+        super(Assembly,self).__init__()
+
+    def execute(self,_obj):
+        self.constraints = None
+        self.buildShape()
+        return False # return False to call LinkBaseExtension::execute()
+
+    def buildShape(self):
+        obj = self.obj
+        if not obj.BuildShape:
+            obj.Shape.nullify()
+            return
+
+        import Part
+        shape = []
+        partGroup = self.getPartGroup(obj)
+        group = partGroup.Group
+        if not group:
+            raise RuntimeError('no parts')
+        if obj.BuildShape == BuildShapeEnum.Cut:
+            shape = Part.getShape(group[0]).Solids
+            if not shape:
+                raise RuntimeError('First part has no solid')
+            if len(shape)>1:
+                shape = [shape[0].fuse(shape[1:])]
+            group = group[1:]
+
+        for o in group:
+            if obj.isElementVisible(o.Name):
+                shape += Part.getShape(o).Solids
+        if not shape:
+            raise RuntimeError('No solids found in parts')
+        if len(shape) == 1:
+            obj.Shape = shape[0]
+        elif obj.BuildShape == BuildShapeEnum.Fuse:
+            obj.Shape = shape[0].fuse(shape[1:])
+        elif obj.BuildShape == BuildShapeEnum.Cut:
+            if len(shape)>2:
+                obj.Shape = shape[0].cut(shape[1].fuse(shape[2:]))
+            else:
+                obj.Shape = shape[0].cut(shape[1])
+        else:
+            obj.Shape = Part.makeCompound(shape)
+
+    def attach(self, obj):
+        obj.addProperty("App::PropertyEnumeration","BuildShape","Base",'')
+        obj.BuildShape = BuildShapeNames
+        super(Assembly,self).attach(obj)
+
+    def linkSetup(self,obj):
+        obj.configLinkProperty('Placement')
+        super(Assembly,self).linkSetup(obj)
+        self.onChanged(obj,'BuildShape')
+
+        # make sure all children are there, first constraint group, then element
+        # group, and finally part group. Call getPartGroup below will make sure
+        # all groups exist. The order of the group is important to make sure
+        # correct rendering and picking behavior
+        self.getPartGroup(True)
+
+    def onChanged(self, obj, prop):
+        if prop == 'BuildShape':
+            if not obj.BuildShape or obj.BuildShape == BuildShapeEnum.Compound:
+                obj.setPropertyStatus('Shape','-Transient')
+            else:
+                obj.setPropertyStatus('Shape','Transient')
+
+    def getConstraintGroup(self, create=False):
+        obj = self.obj
+        try:
+            ret = obj.Group[0]
+            checkType(ret,AsmConstraintGroup)
+            parent = getattr(ret.Proxy,'parent',None)
+            if not parent:
+                ret.Proxy.parent = self
+            elif parent!=self:
+                raise RuntimeError('invalid parent of constraint group '
+                    '{}'.format(objName(ret)))
+            return ret
+        except IndexError:
+            if not create:
+                return # constraint group is optional, so, no exception
+            if obj.Group:
+                raise RuntimeError('Invalid assembly')
+            ret = AsmConstraintGroup.make(obj)
+            obj.setLink({0:ret})
+            return ret
+
+    def getConstraints(self,refresh=False):
+        if not refresh:
+            ret = getattr(self,'constraints',None)
+            if ret:
+                return ret
+        self.constraints = None
+        cstrGroup = self.getConstraintGroup()
+        if not cstrGroup:
+            return
+        ret = []
+        for o in cstrGroup.Group:
+            checkType(o,AsmConstraint)
+            if not o._Type:
+                logger.debug('skip constraint "{}" type '
+                    '"{}"'.format(objName(o),o.Type))
+                continue
+            ret.append(o)
+        self.constraints = ret
+        return self.constraints
+
+    def getElementGroup(self,create=False):
+        obj = self.obj
+        try:
+            ret = obj.Group[1]
+            checkType(ret,AsmElementGroup)
+            parent = getattr(ret.Proxy,'parent',None)
+            if not parent:
+                ret.Proxy.parent = self
+            elif parent!=self:
+                raise RuntimeError('invalid parent of element group '
+                    '{}'.format(objName(ret)))
+            return ret
+        except IndexError:
+            if not create:
+                raise RuntimeError('Missing element group')
+            self.getConstraintGroup(True)
+            ret = AsmElementGroup.make(obj)
+            obj.setLink({1:ret})
+            return ret
+
+    def getPartGroup(self,create=False):
+        obj = self.obj
+        try:
+            ret = obj.Group[2]
+            checkType(ret,AsmPartGroup)
+            parent = getattr(ret.Proxy,'parent',None)
+            if not parent:
+                ret.Proxy.parent = self
+            elif parent!=self:
+                raise RuntimeError(
+                        'invalid parent of part group {}'.format(objName(ret)))
+            return ret
+        except IndexError:
+            if not create:
+                raise RuntimeError('Missing part group')
+            self.getConstraintGroup(True)
+            self.getElementGroup(True)
+            ret = AsmPartGroup.make(obj)
+            obj.setLink({2:ret})
+            return ret
+
+    @staticmethod
+    def make(doc=None,name='Assembly'):
+        if not doc:
+            doc = FreeCAD.ActiveDocument
+        obj = doc.addObject(
+                "Part::FeaturePython",name,Assembly(),None,True)
+        ViewProviderAssembly(obj.ViewObject)
+        obj.Visibility = True
+        return obj
+
+    Info = namedtuple('AssemblyInfo',('Assembly','Object','Subname'))
+
+    @staticmethod
+    def findChild(obj,subname,childType=None,
+            recursive=False,relativeToChild=True):
+        '''
+        Find the immediate child of the first Assembly referenced in 'subs'
+
+        obj: the parent object
+
+        subname: '.' separted sub-object reference, or string list of sub-object
+                 names. Must contain no sub element name.
+
+        childType: optional checking of the child type.
+
+        recursive: If True, continue finding the child of the next assembly.
+
+        relativeToChild: If True, the returned subname is realtive to the child
+        object found, or else, it is relative to the assembly, i.e., including
+        the child's name
+
+        Return None if not found, or (assembly,child,sub), where 'sub' is the
+        remaining sub name list. If recursive is True, then return a list of
+        tuples
+        '''
+        assembly = None
+        child = None
+        idx = -1
+        if isTypeOf(obj,Assembly,True):
+            assembly = obj
+        subs = subname if isinstance(subname,list) else subname.split('.')
+        for i,name in enumerate(subs[:-2]):
+            obj = obj.getSubObject(name+'.',1)
+            if not obj:
+                raise RuntimeError('Cannot find sub object {}'.format(name))
+            if assembly and isTypeOf(obj,childType):
+                child = obj
+                if relativeToChild:
+                    idx = i+1
+                else:
+                    idx = i
+                break
+            assembly = obj if isTypeOf(obj,Assembly,True) else None
+
+        if not child:
+            return
+
+        subs = subs[idx:]
+        ret = Assembly.Info(assembly,child,'.'.join(subs))
+        if not recursive:
+            return ret
+
+        nret = Assembly.findChild(child,subs,childType,True)
+        if nret:
+            return [ret] + nret
+        return [ret]
+
+    @staticmethod
+    def findPartGroup(obj,subname='2.',recursive=False,relativeToChild=True):
+        return Assembly.findChild(
+                obj,subname,AsmPartGroup,recursive,relativeToChild)
+
+    @staticmethod
+    def findElementGroup(obj,subname='1.',relativeToChild=True):
+        return Assembly.findChild(
+                obj,subname,AsmElementGroup,False,relativeToChild)
+
+    @staticmethod
+    def findConstraintGroup(obj,subname='0.',relativeToChild=True):
+        return Assembly.findChild(
+                obj,subname,AsmConstraintGroup,False,relativeToChild)
+
+
+class ViewProviderAssembly(ViewProviderAsmGroup):
+
+    def canDragObject(self,_child):
+        return False
+
+    def canDragObjects(self):
+        return False
+
+    def canDropObject(self,_child):
+        return False
+
+    def canDropObjects(self):
+        return False
+
diff --git a/constraint.py b/constraint.py
new file mode 100644
index 0000000..7af10d6
--- /dev/null
+++ b/constraint.py
@@ -0,0 +1,521 @@
+from collections import namedtuple
+import FreeCAD, FreeCADGui
+import asm3.utils as utils
+import asm3.slvs as slvs
+from asm3.utils import logger, objName
+
+Types = []
+TypeMap = {}
+TypeNameMap = {}
+
+class ConstraintType(type):
+    def __init__(cls, name, bases, attrs):
+        super(ConstraintType,cls).__init__(name,bases,attrs)
+        if cls._id >= 0:
+            if cls._id in TypeMap:
+                raise RuntimeError(
+                        'Duplicate constriant type id {}'.format(cls._id))
+            if cls.slvsFunc():
+                TypeMap[cls._id] = cls
+                TypeNameMap[cls.getName()] = cls
+                cls._idx = len(Types)
+                logger.debug('register constraint "{}":{},{}'.format(
+                    cls.getName(),cls._id,cls._idx))
+                Types.append(cls)
+
+
+# PartName: text name of the part
+# Placement: the original placement of the part
+# Params: 7 parameters that defines the transformation
+# Workplane: a tuple of three entity handles, that is the workplane, the origin
+#            point, and the normal. The workplane, defined by the origin and
+#            norml, is essentially the XY reference plane of the part.
+# EntityMap: string -> entity handle map, for caching
+PartInfo = namedtuple('SolverPartInfo', 
+        ('PartName','Placement','Params','Workplane','EntityMap'))
+
+def _addEntity(etype,system,partInfo,key,shape):
+    key += '.{}'.format(etype)
+    h = partInfo.EntityMap.get(key,None)
+    if h:
+        logger.debug('cache {}: {}'.format(key,h))
+        return h
+    if etype == 'p': # point
+        v = utils.getElementPos(shape)
+        e = system.addPoint3dV(*v)
+    elif etype == 'n': # normal
+        v = utils.getElementNormal(shape)
+        e = system.addNormal3dV(*v)
+    else:
+        raise RuntimeError('unknown entity type {}'.format(etype))
+    h = system.addTransform(e,*partInfo.Params)
+    logger.debug('{}: {},{}, {}'.format(key,h,e,v))
+    partInfo.EntityMap[key] = h
+    return h
+
+def _p(system,partInfo,key,shape):
+    'return a slvs handle of a transformed point derived from "shape"'
+    if not system:
+        if utils.hasCenter(shape):
+            return
+        return 'a vertex or circular edge/face'
+    return _addEntity('p',system,partInfo,key,shape)
+
+def _n(system,partInfo,key,shape):
+    'return a slvs handle of a transformed normal derived from "shape"'
+    if not system:
+        if utils.isAxisOfPlane(shape):
+            return
+        return 'an edge or face with a surface normal'
+    return _addEntity('n',system,partInfo,key,shape)
+
+def _l(system,partInfo,key,shape,retAll=False):
+    'return a pair of slvs handle of the end points of an edge in "shape"'
+    if not system:
+        if utils.isLinearEdge(shape):
+            return
+        return 'a linear edge'
+    key += '.l'
+    h = partInfo.EntityMap.get(key,None)
+    if h:
+        logger.debug('cache {}: {}'.format(key,h))
+    else:
+        v = shape.Edges[0].Vertexes
+        p1 = system.addPoint3dV(*v[0].Point)
+        p2 = system.addPoint3dV(*v[-1].Point)
+        h = system.addLine(p1,p2)
+        h = (h,p1,p2)
+        logger.debug('{}: {}'.format(key,h))
+        partInfo.EntityMap[key] = h
+    return h if retAll else h[0]
+
+def _w(system,partInfo,key,shape,retAll=False):
+    'return a slvs handle of a transformed plane/workplane from "shape"'
+    if not system:
+        if utils.isAxisOfPlane(shape):
+            return
+        return 'an edge or face with a planar surface'
+
+    key2 = key+'.w'
+    h = partInfo.EntityMap.get(key2,None)
+    if h:
+        logger.debug('cache {}: {}'.format(key,h))
+    else:
+        p = _p(system,partInfo,key,shape)
+        n = _n(system,partInfo,key,shape)
+        h = system.addWorkplane(p,n)
+        h = (h,p,n)
+        logger.debug('{}: {}'.format(key,h))
+        partInfo.EntityMap[key2] = h
+    return h if retAll else h[0]
+
+def _c(system,partInfo,key,shape,requireArc=False):
+    'return a slvs handle of a transformed circle/arc derived from "shape"'
+    if not system:
+        r = utils.getElementCircular(shape)
+        if not r or (requireArc and not isinstance(r,list,tuple)):
+            return
+        return 'an cicular arc edge' if requireArc else 'a circular edge'
+    key2 = key+'.c'
+    h = partInfo.EntityMap.get(key2,None)
+    if h:
+        logger.debug('cache {}: {}'.format(key,h))
+    else:
+        h = _w(system,partInfo,key,shape,True)
+        r = utils.getElementCircular(shape)
+        if not r:
+            raise RuntimeError('shape is not cicular')
+        if isinstance(r,(list,tuple)):
+            l = _l(system,partInfo,key,shape,True)
+            h += l[1:]
+            h = system.addArcOfCircleV(*h)
+        elif requireArc:
+            raise RuntimeError('shape is not an arc')
+        else:
+            h = h[1:]
+            h.append(system.addDistanceV(r))
+            h = system.addCircle(*h)
+        logger.debug('{}: {}, {}'.format(key,h,r))
+        partInfo.EntityMap[key2] = h
+    return h
+
+def _a(system,partInfo,key,shape):
+    return _c(system,partInfo,key,shape,True)
+
+
+_PropertyDistance = ('Value','Distance','PropertyDistance','Constraint')
+_PropertyAngle = ('Value','Angle','PropertyAngle','Constraint')
+_PropertyRatio = (None,'Ratio','PropertyFloat','Constraint')
+_PropertyDifference = (None,'Difference','PropertyFloat','Constraint')
+_PropertyDiameter = (None,'Diameter','PropertyFloat','Constraint')
+_PropertyRadius = (None,'Radius','PropertyFloat','Constraint')
+_PropertySupplement = (None,'Supplement','PropertyBool','Constraint',
+        'If True, then the second angle is calculated as 180-angle')
+_PropertyAtEnd = (None,'AtEnd','PropertyBool','Constraint',
+        'If True, then tangent at the end point, or else at the start point')
+
+_ordinal = ['1st', '2nd', '3rd', '4th', '5th', '6th', '7th' ]
+
+class Base:
+    __metaclass__ = ConstraintType
+
+    _id = -1
+    _entities = []
+    _workplane = False
+    _props = []
+    _func = None
+
+    @classmethod
+    def getName(cls):
+        return cls.__name__
+
+    @classmethod
+    def slvsFunc(cls):
+        try:
+            if not cls._func:
+                cls._func = getattr(slvs.System,'add'+cls.getName())
+            return cls._func
+        except AttributeError:
+            logger.error('Invalid slvs constraint "{}"'.format(cls.getName()))
+
+    @classmethod
+    def getEntityDef(cls,group,checkCount,name=None):
+        entities = cls._entities
+        if len(group) != len(entities):
+            if not checkCount and len(group)<len(entities):
+                return entities[:len(group)]
+            if cls._workplane and len(group)==len(entities)+1:
+                entities = list(entities)
+                entities.append(_w)
+            else:
+                if not name:
+                    name = cls.getName()
+                else:
+                    name += ' of type "{}"'.format(cls.getName)
+                raise RuntimeError('Constraint {} has wrong number of '
+                    'elements {}, expecting {}'.format(
+                        name,len(group),len(entities)))
+        return entities
+
+    @classmethod
+    def check(cls,group):
+        entities = cls.getEntityDef(group,False)
+        for i,e in enumerate(entities):
+            o = group[i]
+            msg = e(None,None,None,o)
+            if not msg:
+                continue
+            if i == len(cls._entities):
+                raise RuntimeError('Constraint {} requires the optional {} '
+                    'element to be a planar face for defining a '
+                    'workplane'.format(cls.getName(), _ordinal[i], msg))
+            raise RuntimeError('Constraint {} requires the {} element to be'
+                    ' {}'.format(cls.getName(), _ordinal[i], msg))
+
+    def __init__(self,obj,_props):
+        if obj._Type != self._id:
+            if self._id < 0:
+                raise RuntimeError('invalid constraint type {} id: '
+                    '{}'.format(self.__class__,self._id))
+            obj._Type = self._id
+        for prop in self.__class__._props:
+            obj.addProperty(*prop[1:])
+
+    @classmethod
+    def detach(cls,obj):
+        for prop in cls._props:
+            obj.removeProperty(prop[1])
+
+    def onChanged(self,obj,prop):
+        pass
+
+    @classmethod
+    def getEntities(cls,obj,solver):
+        '''maps fcad element shape to slvs entities'''
+        ret = []
+        for prop in cls._props:
+            v = getattr(obj,prop[1])
+            if prop[0]:
+                v = getattr(v,prop[0])()
+            ret.append(v)
+
+        elements = obj.Proxy.getElements()
+        entities = cls.getEntityDef(elements,True,objName(obj))
+        ret = []
+        for e,o in zip(entities,elements):
+            info = o.Proxy.getInfo()
+            partInfo = solver.getPartInfo(info)
+            ret.append(e(solver.system,partInfo,info.Subname,info.Shape))
+        logger.debug('{}: {}, {}'.format(objName(obj),obj.Type,ret))
+        return ret
+
+    @classmethod
+    def prepare(cls,obj,solver):
+        e = cls.getEntities(obj,solver)
+        cls._func(solver.system,*e,group=solver.group)
+
+
+class Disabled(Base):
+    _id = 0
+    _func = True
+
+    @classmethod
+    def prepare(cls,_obj,_solver):
+        pass
+
+
+class PointsCoincident(Base):
+    _id = 1
+    _entities = (_p,_p)
+    _workplane = True
+
+
+class SameOrientation(Base):
+    _id = 2
+    _entities = (_n,_n)
+
+
+class PointInPlane(Base):
+    _id = 3
+    _entities = (_p,_w)
+
+
+class PointOnLine(Base):
+    _id = 4
+    _entities = (_p,_l)
+    _workplane = True
+
+
+class PointsDistance(Base):
+    _id = 5
+    _entities = (_p,_p)
+    _workplane = True
+    _props = [_PropertyDistance]
+
+
+class PointsProjectDistance(Base):
+    _id = 6
+    _entities = (_p,_p,_l)
+    _props = [_PropertyDistance]
+
+
+class PointPlaneDistance(Base):
+    _id = 7
+    _entities = (_p,_w)
+    _props = [_PropertyDistance]
+
+
+class PointLineDistance(Base):
+    _id = 8
+    _entities = (_p,_l)
+    _workplane = True
+    _props = [_PropertyDistance]
+
+
+class EqualLength(Base):
+    _id = 9
+    _entities = (_l,_l)
+    _workplane = True
+
+
+class LengthRatio(Base):
+    _id = 10
+    _entities = (_l,_l)
+    _workplane = True
+    _props = [_PropertyRatio]
+
+
+class LengthDifference(Base):
+    _id = 11
+    _entities = (_l,_l)
+    _workplane = True
+    _props = [_PropertyDifference]
+
+
+class EqualLengthPointLineDistance(Base):
+    _id = 12
+    _entities = (_p,_l,_l)
+    _workplane = True
+
+
+class EqualPointLineDistance(Base):
+    _id = 13
+    _entities = (_p,_l,_p,_l)
+    _workplane = True
+
+
+class EqualAngle(Base):
+    _id = 14
+    _entities = (_l,_l,_l,_l)
+    _workplane = True
+    _props = [_PropertySupplement]
+
+
+class EqualLineArcLength(Base):
+    _id = 15
+    _entities = (_l,_a)
+    _workplane = True
+
+
+class Symmetric(Base):
+    _id = 16
+    _entities = (_p,_p,_w)
+    _workplane = True
+
+
+class SymmetricHorizontal(Base):
+    _id = 17
+    _entities = (_p,_p,_w)
+
+
+class SymmetricVertical(Base):
+    _id = 18
+    _entities = (_p,_p,_w)
+
+
+class SymmetricLine(Base):
+    _id = 19
+    _entities = (_p,_p,_l,_w)
+
+
+class MidPoint(Base):
+    _id = 20
+    _entities = (_p,_p,_l)
+    _workplane = True
+
+
+class PointsHorizontal(Base):
+    _id = 21
+    _entities = (_p,_p)
+    _workplane = True
+
+
+class PointsVertical(Base):
+    _id = 22
+    _entities = (_p,_p)
+    _workplane = True
+
+
+class LineHorizontal(Base):
+    _id = 23
+    _entities = [_l]
+    _workplane = True
+
+
+class LineVertical(Base):
+    _id = 24
+    _entities = [_l]
+    _workplane = True
+
+
+class Diameter(Base):
+    _id = 25
+    _entities = [_c]
+    _prop = [_PropertyDiameter]
+
+
+class PointOnCircle(Base):
+    _id = 26
+    _entities = [_p,_c]
+
+
+class Angle(Base):
+    _id = 27
+    _entities = (_l,_l)
+    _workplane = True
+    _props = [_PropertyAngle,_PropertySupplement]
+
+
+class Perpendicular(Base):
+    _id = 28
+    _entities = (_l,_l)
+    _workplane = True
+
+
+class Parallel(Base):
+    _id = 29
+    _entities = (_l,_l)
+    _workplane = True
+
+
+class ArcLineTangent(Base):
+    _id = 30
+    _entities = (_c,_l)
+    _props = [_PropertyAtEnd]
+
+
+#  class CubicLineTangent(Base):
+#      _id = 31
+#
+#
+#  class CurvesTangent(Base):
+#      _id = 32
+
+
+class EqualRadius(Base):
+    _id = 33
+    _entities = (_c,_c)
+    _props = [_PropertyRadius]
+
+
+class WhereDragged(Base):
+    _id = 34
+    _entities = [_p]
+    _workplane = True
+
+
+TypeEnum = namedtuple('AsmConstraintEnum',
+        (c.getName() for c in Types))(*range(len(Types)))
+
+def attach(obj,checkType=True):
+    props = None
+    if checkType:
+        props = obj.PropertiesList
+        if not '_Type' in props:
+            raise RuntimeError('Object "{}" has no _Type property'.format(
+                objName(obj)))
+        if 'Type' in props:
+            raise RuntimeError('Object {} already as property "Type"'.format(
+                objName(obj)))
+
+        # The 'Type' property here is to let user select the type in property
+        # editor. It is marked as 'transient' to avoid having to save the
+        # enumeration value for each object.
+        obj.addProperty("App::PropertyEnumeration","Type","Constraint",'',2)
+        obj.Type = TypeEnum._fields
+        idx = 0
+        try:
+            idx = TypeMap[obj._Type]._idx
+        except AttributeError:
+            logger.warn('{} has unknown constraint type {}'.format(
+                objName(obj),obj._Type))
+        obj.Type = idx
+
+    constraintType = TypeNameMap[obj.Type]
+    cstr = getattr(obj.Proxy,'_cstr',None)
+    if type(cstr) is not constraintType:
+        if cstr:
+            cstr.detach(obj)
+        if not props:
+            props = obj.PropertiesList
+        obj.Proxy._cstr = constraintType(obj,props)
+
+
+def onChanged(obj,prop):
+    if prop == 'Type':
+        attach(obj,False)
+        return
+    elif prop == '_Type':
+        obj.Type = TypeMap[obj._Type]._idx
+        return
+    cstr = getattr(obj.Proxy,'_cstr',None)
+    if cstr:
+        cstr.onChanged(obj,prop)
+
+
+def check(tp,group):
+    TypeMap[tp].check(group)
+
+def prepare(cstr,solver):
+    cstr.Proxy._cstr.prepare(cstr,solver)
+
diff --git a/solver.py b/solver.py
new file mode 100644
index 0000000..5184959
--- /dev/null
+++ b/solver.py
@@ -0,0 +1,153 @@
+import FreeCAD, FreeCADGui
+import asm3.slvs as slvs
+import asm3.assembly as asm
+from asm3.utils import logger, objName, isSamePlacement
+import asm3.constraint as constraint
+
+class AsmSolver(object):
+    def __init__(self,assembly,reportFailed):
+        cstrs = assembly.Proxy.getConstraints()
+        if not cstrs:
+            logger.debug('no constraint found in assembly '
+                '{}'.format(objName(assembly)))
+            return
+
+        parts = assembly.Proxy.getPartGroup().Group
+        if len(parts)<=1:
+            logger.debug('not enough parts in {}'.format(objName(assembly)))
+            return
+
+        self.system = slvs.System()
+        self._fixedGroup = 2
+        self.system.GroupHandle = 3 # the import group
+        self.group = 1 # the solving group
+        self._partMap = {}
+        self._cstrMap = {}
+        self._entityMap = {}
+        self._fixedPart = parts[0]
+
+        for cstr in cstrs:
+            logger.debug('preparing {}, type {}'.format(
+                objName(cstr),cstr.Type))
+            self._cstrMap[constraint.prepare(cstr,self)] = cstr
+
+        logger.debug('solving {}'.format(objName(assembly)))
+        ret = self.system.solve(group=self.group,reportFailed=reportFailed)
+        if ret:
+            if reportFailed:
+                msg = 'List of failed constraint:'
+                for f in self.system.Failed:
+                    msg += '\n' + objName(self._cstrMap[f])
+                logger.error(msg)
+            raise RuntimeError('Failed to solve the constraints ({}) '
+                    'in {}'.format(ret,objName(assembly)))
+
+        logger.debug('done sloving, dof {}'.format(self.system.Dof))
+
+        undoDocs = set()
+        for part,partInfo in self._partMap.items():
+            if part == self._fixedPart:
+                continue
+            params = [ self.system.getParam(h).val for h in partInfo.Params ]
+            p = params[:3]
+            q = [params[4],params[5],params[6],params[3]]
+            pla = FreeCAD.Placement(FreeCAD.Vector(*p),FreeCAD.Rotation(*q))
+            if isSamePlacement(partInfo.Placement,pla):
+                logger.debug('not moving {}'.format(partInfo.PartName))
+            else:
+                logger.debug('moving {} {} {} {}'.format(
+                    partInfo.PartName,partInfo.Params,params,pla))
+                asm.AsmElementLink.setPlacement(part,pla,undoDocs)
+
+        for doc in undoDocs:
+            doc.commitTransaction()
+
+
+    def getPartInfo(self,info):
+        partInfo = self._partMap.get(info.Part,None)
+        if partInfo:
+            return partInfo
+
+        # info.Object below is supposed to be the actual part geometry, while
+        # info.Part may be a link to that object. We use info.Object as a key so
+        # that multiple info.Part can share the same entity map.
+        #
+        # TODO: It is actually more complicated than that. Becuase info.Object
+        # itself may be a link, and followed by a chain of multiple links. It's
+        # complicated because each link can either have a linked placement or
+        # not, depends on its LinkTransform property, meaning that their
+        # Placement may be chained or independent.  Ideally, we should explode
+        # the link chain, and recreate the transformation dependency using slvs
+        # transfomation entity. We'll leave that for another day, maybe... 
+        #
+        # The down side for now is that we may have redundant entities, and
+        # worse, the solver may not be able to get correct result if there are
+        # placement dependent links among parts.  So, one should avoid using
+        # LinkTransform in most cases.
+        #
+        entityMap = self._entityMap.get(info.Object,None)
+        if not entityMap:
+            entityMap = {}
+            self._entityMap[info.Object] = entityMap
+
+        if info.Part == self._fixedPart:
+            g = self._fixedGroup
+        else:
+            g = self.group
+        q = info.Placement.Rotation.Q
+        vals = list(info.Placement.Base) + [q[3],q[0],q[1],q[2]]
+        params = [self.system.addParamV(v,g) for v in vals]
+
+        p = self.system.addPoint3d(*params[:3],group=g)
+        n = self.system.addNormal3d(*params[3:],group=g)
+        w = self.system.addWorkplane(p,n,group=g)
+        h = (w,p,n)
+
+        logger.debug('{} {}, {}, {}, {}'.format(
+            info.PartName,info.Placement,h,params,vals))
+
+        partInfo = constraint.PartInfo(
+                info.PartName, info.Placement.copy(),params,h,entityMap)
+        self._partMap[info.Part] = partInfo
+        return partInfo
+
+
+def solve(objs=None,recursive=True,reportFailed=True,recompute=True):
+    if not objs:
+        objs = FreeCAD.ActiveDocument.Objects
+    elif not isinstance(objs,(list,tuple)):
+        objs = [objs]
+
+    if not objs:
+        logger.error('no objects')
+        return
+
+    if recompute:
+        docs = set()
+        for o in objs:
+            docs.add(o.Document)
+        for d in docs:
+            logger.debug('recomputing {}'.format(d.Name))
+            d.recompute()
+
+    if recursive:
+        # Get all dependent object, including external ones, and return as a
+        # topologically sorted list.
+        objs = FreeCAD.getDependentObjects(objs,False,True)
+
+    assemblies = []
+    for obj in objs:
+        if asm.isTypeOf(obj,asm.Assembly):
+            logger.debug('adding assembly {}'.format(objName(obj)))
+            assemblies.append(obj)
+
+    if not assemblies:
+        logger.error('no assembly found')
+        return
+
+    for assembly in assemblies:
+        logger.debug('solving assembly {}'.format(objName(assembly)))
+        AsmSolver(assembly,reportFailed)
+        if recompute:
+            assembly.Document.recompute()
+
diff --git a/utils.py b/utils.py
new file mode 100644
index 0000000..b26aa3f
--- /dev/null
+++ b/utils.py
@@ -0,0 +1,341 @@
+'''
+Collection of helper function to extract geometry properties from OCC elements
+
+Most of the functions are borrowed directly from assembly2lib.py or lib3D.py in
+assembly2
+'''
+
+import FreeCAD, FreeCADGui, Part
+import numpy
+
+import asm3.FCADLogger
+logger = asm3.FCADLogger.FCADLogger('assembly3')
+
+def objName(obj):
+    if obj.Label == obj.Name:
+        return obj.Name
+    return '{}({})'.format(obj.Name,obj.Label)
+
+def isLine(param):
+    if hasattr(Part,"LineSegment"):
+        return isinstance(param,(Part.Line,Part.LineSegment))
+    else:
+        return isinstance(param,Part.Line)
+
+def getElement(obj,tp):
+    if isinstance(obj,tuple):
+       obj = obj[0].getSubObject(obj[1])
+    if not isinstance(obj,Part.Shape):
+        return
+    if obj.isNull() or not tp:
+        return
+    if tp == 'Vertex':
+        vertexes = obj.Vertexes
+        if len(vertexes)==1 and not obj.Edges:
+            return vertexes[0]
+    elif tp == 'Edge':
+        edges = obj.Edges
+        if len(edges)==1 and not obj.Faces:
+            return edges[0]
+    elif tp == 'Face':
+        faces = obj.Faces
+        if len(faces)==1:
+            return faces[0]
+
+def isPlane(obj):
+    face = getElement(obj,'Face')
+    if not face:
+        return False
+    elif str(face.Surface) == '<Plane object>':
+        return True
+    elif hasattr(face.Surface,'Radius'):
+        return False
+    elif str(face.Surface).startswith('<SurfaceOfRevolution'):
+        return False
+    else:
+        _plane_norm,_plane_pos,error = fit_plane_to_surface1(face.Surface)
+        error_normalized = error / face.BoundBox.DiagonalLength
+        return error_normalized < 10**-6
+
+def isCylindricalPlane(obj):
+    face = getElement(obj,'Face')
+    if not face:
+        return False
+    elif hasattr(face.Surface,'Radius'):
+        return True
+    elif str(face.Surface).startswith('<SurfaceOfRevolution'):
+        return True
+    elif str(face.Surface) == '<Plane object>':
+        return False
+    else:
+        _axis,_center,error=fit_rotation_axis_to_surface1(face.Surface)
+        error_normalized = error / face.BoundBox.DiagonalLength
+        return error_normalized < 10**-6
+
+def isAxisOfPlane(obj):
+    face = getElement(obj,'Face')
+    if not face:
+        return False
+    if str(face.Surface) == '<Plane object>':
+        return True
+    else:
+        _axis,_center,error=fit_rotation_axis_to_surface1(face.Surface)
+        error_normalized = error / face.BoundBox.DiagonalLength
+        return error_normalized < 10**-6
+
+def isCircularEdge(obj):
+    edge = getElement(obj,'Edge')
+    if not edge:
+        return False
+    elif not hasattr(edge, 'Curve'): #issue 39
+        return False
+    if hasattr( edge.Curve, 'Radius' ):
+        return True
+    elif isLine(edge.Curve):
+        return False
+    else:
+        BSpline = edge.Curve.toBSpline()
+        try:
+            arcs = BSpline.toBiArcs(10**-6)
+        except Exception:  #FreeCAD exception thrown ()
+            return False
+        if all( hasattr(a,'Center') for a in arcs ):
+            centers = numpy.array([a.Center for a in arcs])
+            sigma = numpy.std( centers, axis=0 )
+            return max(sigma) < 10**-6
+        return False
+
+def isLinearEdge(obj):
+    edge = getElement(obj,'Edge')
+    if not edge:
+        return False
+    elif not hasattr(edge, 'Curve'): #issue 39
+        return False
+    if isLine(edge.Curve):
+        return True
+    elif hasattr( edge.Curve, 'Radius' ):
+        return False
+    else:
+        BSpline = edge.Curve.toBSpline()
+        try:
+            arcs = BSpline.toBiArcs(10**-6)
+        except Exception:  #FreeCAD exception thrown ()
+            return False
+        if all(isLine(a) for a in arcs):
+            lines = arcs
+            D = numpy.array([L.tangent(0)[0] for L in lines]) #D(irections)
+            return numpy.std( D, axis=0 ).max() < 10**-9
+        return False
+
+def isVertex(obj):
+    return getElement(obj,'Vertex') is not None
+
+def hasCenter(obj):
+    return isVertex(obj) or isCircularEdge(obj) or \
+            isAxisOfPlane(obj) or isSphericalSurface(obj)
+
+def isSphericalSurface(obj):
+    face = getElement(obj,'Face')
+    if not face:
+        return False
+    return str( face.Surface ).startswith('Sphere ')
+
+def getElementPos(obj):
+    pos = None
+    vertex = getElement(obj,'Vertex')
+    if vertex:
+        return vertex.Point
+    face = getElement(obj,'Face')
+    if face:
+        surface = face.Surface
+        if str(surface) == '<Plane object>':
+            #  pos = face.BoundBox.Center
+            pos = surface.Position
+        elif all( hasattr(surface,a) for a in ['Axis','Center','Radius'] ):
+            pos = surface.Center
+        elif str(surface).startswith('<SurfaceOfRevolution'):
+            pos = face.Edges1.Curve.Center
+        else: #numerically approximating surface
+            _plane_norm, plane_pos, error = \
+                    fit_plane_to_surface1(face.Surface)
+            error_normalized = error / face.BoundBox.DiagonalLength
+            if error_normalized < 10**-6: #then good plane fit
+                pos = plane_pos
+            _axis, center, error = \
+                    fit_rotation_axis_to_surface1(face.Surface)
+            error_normalized = error / face.BoundBox.DiagonalLength
+            if error_normalized < 10**-6: #then good rotation_axis fix
+                pos = center
+    else:
+        edge = getElement(obj,'Edge')
+        if edge:
+            if isLine(edge.Curve):
+                pos = edge.Vertexes[-1].Point
+            elif hasattr( edge.Curve, 'Center'): #circular curve
+                pos = edge.Curve.Center
+            else:
+                BSpline = edge.Curve.toBSpline()
+                arcs = BSpline.toBiArcs(10**-6)
+                if all( hasattr(a,'Center') for a in arcs ):
+                    centers = numpy.array([a.Center for a in arcs])
+                    sigma = numpy.std( centers, axis=0 )
+                    if max(sigma) < 10**-6: #then circular curce
+                        pos = centers[0]
+                elif all(isLine(a) for a in arcs):
+                    lines = arcs
+                    D = numpy.array(
+                            [L.tangent(0)[0] for L in lines]) #D(irections)
+                    if numpy.std( D, axis=0 ).max() < 10**-9: #then linear curve
+                        return lines[0].value(0)
+    return pos
+
+
+def getElementAxis(obj):
+    axis = None
+    face = getElement(obj,'Face')
+    if face:
+        surface = face.Surface
+        if hasattr(surface,'Axis'):
+            axis = surface.Axis
+        elif str(surface).startswith('<SurfaceOfRevolution'):
+            axis = face.Edges[0].Curve.Axis
+        else: #numerically approximating surface
+            plane_norm, _plane_pos, error = \
+                    fit_plane_to_surface1(face.Surface)
+            error_normalized = error / face.BoundBox.DiagonalLength
+            if error_normalized < 10**-6: #then good plane fit
+                axis = plane_norm
+            axis_fitted, _center, error = \
+                    fit_rotation_axis_to_surface1(face.Surface)
+            error_normalized = error / face.BoundBox.DiagonalLength
+            if error_normalized < 10**-6: #then good rotation_axis fix
+                axis = axis_fitted
+    else:
+        edge = getElement(obj,'Edge')
+        if edge:
+            if isLine(edge.Curve):
+                axis = edge.Curve.tangent(0)[0]
+            elif hasattr( edge.Curve, 'Axis'): #circular curve
+                axis =  edge.Curve.Axis
+            else:
+                BSpline = edge.Curve.toBSpline()
+                arcs = BSpline.toBiArcs(10**-6)
+                if all( hasattr(a,'Center') for a in arcs ):
+                    centers = numpy.array([a.Center for a in arcs])
+                    sigma = numpy.std( centers, axis=0 )
+                    if max(sigma) < 10**-6: #then circular curce
+                        axis = arcs[0].Axis
+                if all(isLine(a) for a in arcs):
+                    lines = arcs
+                    D = numpy.array(
+                            [L.tangent(0)[0] for L in lines]) #D(irections)
+                    if numpy.std( D, axis=0 ).max() < 10**-9: #then linear curve
+                        return D[0]
+    return axis
+
+def axisToNormal(v):
+    return FreeCAD.Rotation(FreeCAD.Vector(0,0,1),v).Q
+
+def getElementNormal(obj):
+    return axisToNormal(getElementAxis(obj))
+
+def getElementCircular(obj):
+    'return radius if it is closed, or a list of two endpoints'
+    edge = getElement(obj,'Edge')
+    if not edge:
+        return
+    elif not hasattr(edge, 'Curve'): #issue 39
+        return
+    c = edge.Curve
+    if hasattr( c, 'Radius' ):
+        if edge.Closed:
+            return c.Radius
+    elif isLine(edge.Curve):
+        return
+    else:
+        BSpline = edge.Curve.toBSpline()
+        try:
+            arc = BSpline.toBiArcs(10**-6)[0]
+        except Exception:  #FreeCAD exception thrown ()
+            return
+        if edge.Closed:
+            return arc[0].Radius
+    return [v.Point for v in edge.Vertexes]
+
+def fit_plane_to_surface1( surface, n_u=3, n_v=3 ):
+    'borrowed from assembly2 lib3D.py'
+    uv = sum( [ [ (u,v) for u in numpy.linspace(0,1,n_u)]
+        for v in numpy.linspace(0,1,n_v) ], [] )
+    # positions at u,v points
+    P = [ surface.value(u,v) for u,v in uv ]
+    N = [ numpy.cross( *surface.tangent(u,v) ) for u,v in uv ]
+    # plane's normal, averaging done to reduce error
+    plane_norm = sum(N) / len(N)
+    plane_pos = P[0]
+    error = sum([ abs( numpy.dot(p - plane_pos, plane_norm) ) for p in P ])
+    return plane_norm, plane_pos, error
+
+def fit_rotation_axis_to_surface1( surface, n_u=3, n_v=3 ):
+    '''
+    should work for cylinders and pssibly cones (depending on the u,v mapping)
+
+    borrowed from assembly2 lib3D.py
+    '''
+    uv = sum( [ [ (u,v) for u in numpy.linspace(0,1,n_u)]
+        for v in numpy.linspace(0,1,n_v) ], [] )
+    # positions at u,v points
+    P = [ numpy.array(surface.value(u,v)) for u,v in uv ]
+    N = [ numpy.cross( *surface.tangent(u,v) ) for u,v in uv ]
+    intersections = []
+    for i in range(len(N)-1):
+        for j in range(i+1,len(N)):
+            # based on the distance_between_axes( p1, u1, p2, u2) function,
+            if 1 - abs(numpy.dot( N[i], N[j])) < 10**-6:
+                continue #ignore parrallel case
+            p1_x, p1_y, p1_z = P[i]
+            u1_x, u1_y, u1_z = N[i]
+            p2_x, p2_y, p2_z = P[j]
+            u2_x, u2_y, u2_z = N[j]
+            t1_t1_coef = u1_x**2 + u1_y**2 + u1_z**2 #should equal 1
+            # collect( expand(d_sqrd), [t1*t2] )
+            t1_t2_coef = -2*u1_x*u2_x - 2*u1_y*u2_y - 2*u1_z*u2_z
+            t2_t2_coef = u2_x**2 + u2_y**2 + u2_z**2 #should equal 1 too
+            t1_coef    = 2*p1_x*u1_x + 2*p1_y*u1_y + 2*p1_z*u1_z - \
+                    2*p2_x*u1_x - 2*p2_y*u1_y - 2*p2_z*u1_z
+            t2_coef    =-2*p1_x*u2_x - 2*p1_y*u2_y - 2*p1_z*u2_z + \
+                    2*p2_x*u2_x + 2*p2_y*u2_y + 2*p2_z*u2_z
+            A = numpy.array([ [ 2*t1_t1_coef , t1_t2_coef ],
+                [ t1_t2_coef, 2*t2_t2_coef ] ])
+            b = numpy.array([ t1_coef, t2_coef])
+            try:
+                t1, t2 = numpy.linalg.solve(A,-b)
+            except numpy.linalg.LinAlgError:
+                continue
+            pos_t1 = P[i] + numpy.array(N[i])*t1
+            pos_t2 = P[j] + N[j]*t2
+            intersections.append( pos_t1 )
+            intersections.append( pos_t2 )
+    if len(intersections) < 2:
+        error = numpy.inf
+        return 0, 0, error
+    else: 
+        # fit vector to intersection points; 
+        # http://mathforum.org/library/drmath/view/69103.html
+        X = numpy.array(intersections)
+        centroid = numpy.mean(X,axis=0)
+        M = numpy.array([i - centroid for i in intersections ])
+        A = numpy.dot(M.transpose(), M)
+        # numpy docs: s : (..., K) The singular values for every matrix, 
+        # sorted in descending order.
+        _U,s,V = numpy.linalg.svd(A)
+        axis_pos = centroid
+        axis_dir = V[0]
+        error = s[1] #dont know if this will work
+        return axis_dir, axis_pos, error
+
+_tol = 10e-7
+
+def isSamePlacement(pla1,pla2):
+    return pla1.Base.distanceToPoint(pla2.Base) < _tol and \
+        numpy.norm(numpy.array(pla1.Rotation.Q) - \
+                   numpy.array(pla2.Rotation.Q)) < _tol