From 8a03bc40b38cd665fb59ea2c1d4afddc39af008e Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Mon, 22 Jun 2015 22:44:28 +0200
Subject: [PATCH 1/9] Update feamodel.py

---
 pycalculix/feamodel.py | 68 +++++++++++++++++++++++-------------------
 1 file changed, 38 insertions(+), 30 deletions(-)

diff --git a/pycalculix/feamodel.py b/pycalculix/feamodel.py
index ef0966c..74de570 100644
--- a/pycalculix/feamodel.py
+++ b/pycalculix/feamodel.py
@@ -158,6 +158,22 @@ def set_ediv(self, items, ediv):
         items = self.get_items(items)
         for line in items:
             line.set_ediv(ediv)
+            
+    def set_esize(self, items, esize):
+        """Sets the element size on the passed line.
+
+        Args:
+            items (str or SignLine or SignArc or list): lines or points to set esize on
+
+                - str: 'L0'
+                - list of str ['L0', 'L1', 'P3']
+                - list of SignLine or SignArc part.bottom or part.hole[-1]
+
+            esize (float): size of the mesh elements on the line
+        """
+        items = self.get_items(items)
+        for item in items:
+            item.set_esize(esize)
 
     def set_units(self, dist_unit='m', cfswitch=False):
         """Sets the units that will be displayed when plotting.
@@ -1347,7 +1363,7 @@ def __read_inp(self, fname):
                     L = line.split(',')
                     L = [int(a.strip()) for a in L]
                     enum = L[0]
-                    nlist = [N.idget(a) for a in L[1:]]
+                    nlist = [Dict_NodeIDs[a] for a in L[1:]]
                     e = mesh.Element(enum, etype, nlist)
                     faces = e.faces
                     E.append(e)
@@ -1371,6 +1387,7 @@ def __read_inp(self, fname):
             # mode setting
             if '*Node' in line or '*NODE' in line:
                 mode = 'nmake'
+                print('mode=nmake')
             elif '*Element' in line or '*ELEMENT' in line:
                 L = line.split(',') # split it based on commas
                 e = L[1].split('=')
@@ -1383,6 +1400,7 @@ def __read_inp(self, fname):
                     set_type = 'E'
                     sets[set_type][set_name] = []
                     mode = 'emake'
+                    print('mode=emake')
                 else:
                     mode = None
             elif '*ELSET' in line:
@@ -1392,6 +1410,7 @@ def __read_inp(self, fname):
                 set_type = 'E'
                 sets[set_type][set_name] = []
                 mode = 'set'
+                print('mode=set')
             elif '*NSET' in line:
                 L = line.split(',')
                 e = L[1].split('=')
@@ -1399,8 +1418,11 @@ def __read_inp(self, fname):
                 set_type = 'N'
                 sets[set_type][set_name] = []
                 mode = 'set'
+                print('mode=set')
         f.close()
 
+        print('inputfile read')
+		
         # loop through sets and remove empty sets
         # store sets to delete
         todel = []
@@ -1412,7 +1434,7 @@ def __read_inp(self, fname):
         for adict in todel:
             (set_type, set_name) = (adict['set_type'], adict['set_name'])
             del sets[set_type][set_name]
-            #print('Empty set type:%s name:%s deleted' % (set_type, set_name))
+            print('Empty set type:%s name:%s deleted' % (set_type, set_name))
 
         # this resets the min element to number 1
         if E.get_minid() > 1:
@@ -1492,11 +1514,11 @@ def __read_inp(self, fname):
         print('Nodes: %i' % len(N))
         print('Done reading Calculix/Abaqus .inp file')
 
-    def mesh(self, fineness=1.0, mesher='gmsh'):
+    def mesh(self, element_size=1.0, mesher='gmsh'):
         """Meshes all parts.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
+            element_size (float): 0.0001 - 1.0, how fine the mesh is.
 
                 - Low numbers are very fine, higher numbers are coarser.
             mesher (str): the mesher to use
@@ -1505,15 +1527,15 @@ def mesh(self, fineness=1.0, mesher='gmsh'):
                 - 'cgx': mesh with Calculix cgx, it doesn't allow holes
         """
         if mesher == 'gmsh':
-            self.__mesh_gmsh(fineness)
+            self.__mesh_gmsh(element_size)
         elif mesher == 'cgx':
-            self.__mesh_cgx(fineness)
+            self.__mesh_cgx(element_size)
 
-    def __mesh_gmsh(self, fineness):
+    def __mesh_gmsh(self, element_size):
         """Meshes all parts using the Gmsh mesher.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
+            element_size (float): 0.0001 - 1.0, how fine the mesh is.
 
                 Low numbers are very fine, higher numbers are coarser.
         """
@@ -1524,7 +1546,10 @@ def __mesh_gmsh(self, fineness):
 
         # write all points
         for pt in self.points:
-            txtline = 'Point(%i) = {%f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0)
+            if pt.esize == None:
+                txtline = 'Point(%i) = {%f, %f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0, element_size if self.eshape=='tri' else element_size*2.)
+            else:
+                txtline = 'Point(%i) = {%f, %f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0, pt.esize if self.eshape=='tri' else pt.esize*2.)
             geo.append(txtline)
 
         # start storing an index number
@@ -1544,21 +1569,6 @@ def __mesh_gmsh(self, fineness):
                 txtline = 'Line(%i) = {%i,%i};' % (ind, pt1, pt2)
             geo.append(txtline)
 
-            # set division if we have it
-            if line.ediv != None:
-                ndiv = line.ediv+1
-                esize = line.length()/line.ediv
-                if self.eshape == 'quad':
-                    ndiv = line.ediv/2+1
-                    esize = esize*2
-                    # this is needed because quad recombine
-                    # splits 1 element into 2
-                txtline = 'Transfinite Line{%i} = %i;' % (ind, ndiv)
-                print('LINE ELEMENT SIZE: %f, MAKES %i ELEMENTS'
-                      % (line.length()/line.ediv, line.ediv))
-                geo.append(txtline)
-                geo.append('Characteristic Length {%i,%i} = %f;'
-                           % (pt1, pt2, esize))
             ind += 1
 
         # write all areas
@@ -1613,8 +1623,6 @@ def __mesh_gmsh(self, fineness):
             geo.append(txtline)
 
         # set the meshing options
-        geo.append('Mesh.CharacteristicLengthFactor = '
-                   +str(fineness)+'; //mesh fineness')
         geo.append('Mesh.RecombinationAlgorithm = 1; //blossom')
 
         if self.eshape == 'quad':
@@ -1663,11 +1671,11 @@ def __mesh_gmsh(self, fineness):
         # read in the calculix mesh
         self.__read_inp(self.fname+'.inp')
 
-    def __mesh_cgx(self, fineness):
+    def __mesh_cgx(self, element_size):
         """Meshes all parts using the Calculix cgx mesher.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
+            element_size (float): 0.0001 - 1.0, how fine the mesh is.
 
                 Low numbers are very fine, higher numbers are coarser.
         """
@@ -1693,8 +1701,8 @@ def __mesh_cgx(self, fineness):
         cgx_elements['quad2plstrain'] = 'qu8e'
         cgx_elements['quad1plstrain'] = 'qu4e'
 
-        num = 1.0/fineness
-        emult = int(round(num)) # this converts fineness to mesh multiplier
+        num = 1.0/element_size
+        emult = int(round(num)) # this converts element_size to mesh multiplier
 
         # write all points
         for point in self.points:

From c6b57f09f4a2db5184c492211c5f32eb5ad09b57 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Mon, 22 Jun 2015 22:45:38 +0200
Subject: [PATCH 2/9] Update geometry.py

---
 pycalculix/geometry.py | 41 +++++++++++++++++++++++++++++++++++++----
 1 file changed, 37 insertions(+), 4 deletions(-)

diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py
index 51e6ed5..8bc20cb 100644
--- a/pycalculix/geometry.py
+++ b/pycalculix/geometry.py
@@ -67,6 +67,7 @@ def __init__(self, x, y, z=0):
         self.lines = []
         self.arc_center = False
         self.on_part = True
+        self.esize = None
         base_classes.Idobj.__init__(self)
 
     def get_name(self):
@@ -258,6 +259,10 @@ def __str__(self):
         """Returns string listing object type, id number, and coordinates"""
         val = 'Point %s, (x, y)=(%.3f, %.3f)' % (self.get_name(), self.x, self.y)
         return val
+    
+    def set_esize(self, size):
+        self.esize = size
+        
 
 
 class Line(base_classes.Idobj):
@@ -351,7 +356,20 @@ def set_ediv(self, ediv):
         Args:
             ediv (int): number of required elements on this line
         """
+        print('Funktion ediv aufgerufen')
         self.ediv = ediv
+        self.pt(0).set_esize(self.length()/ediv)
+        self.pt(1).set_esize(self.length()/ediv)
+        
+        
+    def set_esize(self, esize):
+        """Sets the size of mesh elements on the line when meshing.
+
+        Args:
+            esize (float): size of mesh elements on this line
+        """
+        self.pt(0).set_esize(esize)
+        self.pt(1).set_esize(esize)
 
     def signed_copy(self, sign):
         """Returns a SignLine copy of this Line with the passed sign."""
@@ -675,7 +693,11 @@ def set_pt(self, index, point):
     def set_ediv(self, ediv):
         """Applies the element divisions onto the parent line."""
         self.line.set_ediv(ediv)
-
+        
+    def set_esize(self, esize):
+        """Applies the element size onto the parent line."""
+        self.line.set_esize(esize)
+        
     def set_lineloop(self, lineloop):
         """Sets the parent LineLoop"""
         # this is needed to cascade set ediv up to FEA model and down onto
@@ -869,15 +891,26 @@ def save_to_points(self):
         """This stores this line in the line's child points."""
         for point in self.allpoints:
             point.save_line(self)
-
+        
     def set_ediv(self, ediv):
-        """Sets the number of element divisions on the arc when meshing.
+        """Sets the number of element divisions on the line when meshing.
 
         Args:
-            ediv (int): number of required elements on this arc
+            ediv (int): number of required elements on this line
         """
         self.ediv = ediv
+        self.pt(0).set_esize(self.length()/ediv)
+        self.pt(1).set_esize(self.length()/ediv)
+        
+    def set_esize(self, esize):
+        """Sets the size of mesh elements on the line when meshing.
 
+        Args:
+            esize (float): size of mesh elements on this line
+        """
+        self.pt(0).set_esize(esize)
+        self.pt(1).set_esize(esize)
+        
     def signed_copy(self, sign):
         """Returns a SignArc instance of this Arc with the passed sign."""
         return SignArc(self, sign)

From d827fab767481016e0e1eb73db2a0ee0f98fd075 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Mon, 22 Jun 2015 23:08:31 +0200
Subject: [PATCH 3/9] Fixed meshing problems with gmsh

- meshing with transfinite line lead to bad meshes in the intersection between lines
- line next to transfinite line was not affected by mesh size changes (easy to see in dam.py examplte with fine mesh)
- I changed mesh size definition to 'element size at point'
- added function 'set_esize' to define element size instead of numbers
- meshing should also be more stable with this
- the same size definition lead to quad elements be half the size of tri elements, added simple fix for this in __mesh__gmsh function
---
 pycalculix/feamodel.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/pycalculix/feamodel.py b/pycalculix/feamodel.py
index 74de570..53cee80 100644
--- a/pycalculix/feamodel.py
+++ b/pycalculix/feamodel.py
@@ -1387,7 +1387,6 @@ def __read_inp(self, fname):
             # mode setting
             if '*Node' in line or '*NODE' in line:
                 mode = 'nmake'
-                print('mode=nmake')
             elif '*Element' in line or '*ELEMENT' in line:
                 L = line.split(',') # split it based on commas
                 e = L[1].split('=')
@@ -1400,7 +1399,6 @@ def __read_inp(self, fname):
                     set_type = 'E'
                     sets[set_type][set_name] = []
                     mode = 'emake'
-                    print('mode=emake')
                 else:
                     mode = None
             elif '*ELSET' in line:
@@ -1410,7 +1408,6 @@ def __read_inp(self, fname):
                 set_type = 'E'
                 sets[set_type][set_name] = []
                 mode = 'set'
-                print('mode=set')
             elif '*NSET' in line:
                 L = line.split(',')
                 e = L[1].split('=')
@@ -1418,11 +1415,8 @@ def __read_inp(self, fname):
                 set_type = 'N'
                 sets[set_type][set_name] = []
                 mode = 'set'
-                print('mode=set')
         f.close()
-
-        print('inputfile read')
-		
+        
         # loop through sets and remove empty sets
         # store sets to delete
         todel = []
@@ -1434,7 +1428,7 @@ def __read_inp(self, fname):
         for adict in todel:
             (set_type, set_name) = (adict['set_type'], adict['set_name'])
             del sets[set_type][set_name]
-            print('Empty set type:%s name:%s deleted' % (set_type, set_name))
+            #print('Empty set type:%s name:%s deleted' % (set_type, set_name))
 
         # this resets the min element to number 1
         if E.get_minid() > 1:

From bd1bb506f07d543df0c7689bab342db7c1914501 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Mon, 22 Jun 2015 23:09:45 +0200
Subject: [PATCH 4/9] Fixed meshing problems with gmsh

- meshing with transfinite line lead to bad meshes in the intersection between lines
- line next to transfinite line was not affected by mesh size changes (easy to see in dam.py examplte with fine mesh)
- I changed mesh size definition to 'element size at point'
- added function 'set_esize' to define element size instead of numbers
- meshing should also be more stable with this
- the same size definition lead to quad elements be half the size of tri elements, added simple fix for this in __mesh__gmsh function
---
 pycalculix/geometry.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py
index 8bc20cb..eb53522 100644
--- a/pycalculix/geometry.py
+++ b/pycalculix/geometry.py
@@ -356,7 +356,7 @@ def set_ediv(self, ediv):
         Args:
             ediv (int): number of required elements on this line
         """
-        print('Funktion ediv aufgerufen')
+        
         self.ediv = ediv
         self.pt(0).set_esize(self.length()/ediv)
         self.pt(1).set_esize(self.length()/ediv)
@@ -893,20 +893,20 @@ def save_to_points(self):
             point.save_line(self)
         
     def set_ediv(self, ediv):
-        """Sets the number of element divisions on the line when meshing.
+        """Sets the number of element divisions on the arc when meshing.
 
         Args:
-            ediv (int): number of required elements on this line
+            ediv (int): number of required elements on this arc
         """
         self.ediv = ediv
         self.pt(0).set_esize(self.length()/ediv)
         self.pt(1).set_esize(self.length()/ediv)
         
     def set_esize(self, esize):
-        """Sets the size of mesh elements on the line when meshing.
+        """Sets the size of mesh elements on the arc when meshing.
 
         Args:
-            esize (float): size of mesh elements on this line
+            esize (float): size of mesh elements on this arc
         """
         self.pt(0).set_esize(esize)
         self.pt(1).set_esize(esize)

From 8392dbbf4fd84c0d432b28f38a5f33f48d6bdca3 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Sun, 5 Jul 2015 17:55:49 +0200
Subject: [PATCH 5/9] Add mesh mode functionality  (This functionality only
 works with 'gmsh' mesher !): - able to choose between meshmodes 'fineness'
 and 'esize' - fineness: default setting, like before - esize:  - element size
 on points is set                - mode is set if meshmode 'esize' is chosen
 (mesh(meshmode='esize') or function set_esize is called                -
 function set_esize is working for areas, lines, arcs and points (esize is
 cascaded to size on points)                - if the elementsize is set
 multiple times, the prior defined esize properties gets overwritten          
      - if mode esize is chosen but ediv is set to lines, ediv is transformed
 to esize (and prior defined esize gets overwritten, therefore mixing those is
 not recommended)

---
 pycalculix/feamodel.py | 128 ++++++++++++++++++++++++++++++++++-------
 pycalculix/geometry.py |  48 +++++++++++++++-
 2 files changed, 152 insertions(+), 24 deletions(-)

diff --git a/pycalculix/feamodel.py b/pycalculix/feamodel.py
index ef0966c..6bd38fb 100644
--- a/pycalculix/feamodel.py
+++ b/pycalculix/feamodel.py
@@ -158,6 +158,22 @@ def set_ediv(self, items, ediv):
         items = self.get_items(items)
         for line in items:
             line.set_ediv(ediv)
+            
+    def set_esize(self, items, esize):
+        """Sets the element size on the passed line.
+
+        Args:
+            items (str or SignLine or SignArc or list): lines or points to set esize on
+
+                - str: 'L0'
+                - list of str ['L0', 'L1', 'P3']
+                - list of SignLine or SignArc part.bottom or part.hole[-1]
+
+            esize (float): size of the mesh elements on the line
+        """
+        items = self.get_items(items)
+        for item in items:
+            item.set_esize(esize)
 
     def set_units(self, dist_unit='m', cfswitch=False):
         """Sets the units that will be displayed when plotting.
@@ -1347,7 +1363,7 @@ def __read_inp(self, fname):
                     L = line.split(',')
                     L = [int(a.strip()) for a in L]
                     enum = L[0]
-                    nlist = [N.idget(a) for a in L[1:]]
+                    nlist = [Dict_NodeIDs[a] for a in L[1:]]
                     e = mesh.Element(enum, etype, nlist)
                     faces = e.faces
                     E.append(e)
@@ -1400,7 +1416,7 @@ def __read_inp(self, fname):
                 sets[set_type][set_name] = []
                 mode = 'set'
         f.close()
-
+        
         # loop through sets and remove empty sets
         # store sets to delete
         todel = []
@@ -1492,30 +1508,78 @@ def __read_inp(self, fname):
         print('Nodes: %i' % len(N))
         print('Done reading Calculix/Abaqus .inp file')
 
-    def mesh(self, fineness=1.0, mesher='gmsh'):
+    def mesh(self, size=1.0, meshmode='fineness', mesher='gmsh'):
         """Meshes all parts.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
-
-                - Low numbers are very fine, higher numbers are coarser.
+            size (float): 
+            
+                - if meshmode == 'fineness' (default):
+                    - mesh size is adapted to geometry size
+                    - set size = 0.0001 - 1.0, to define how fine the mesh is.
+                    - Low numbers are very fine, higher numbers are coarser.
+                
+                - if meshmode == 'esize':
+                    - element size is kept constant
+                    - choose it depending on geometry size
+                    - it should be reduced e.g. at arcs with small radius, by calling line.esize function
+                
+            meshmode (str):
+                
+                - 'fineness': adapt mesh size to geometry
+                - 'esize': keep explicitly defined element size
+                
+                meshmode is changed to 'esize' is used if esize property is set to points or lines
+            
             mesher (str): the mesher to use
 
                 - 'gmsh': mesh with Gmsh, this is reccomended, it allows holes
                 - 'cgx': mesh with Calculix cgx, it doesn't allow holes
         """
+        
+        
+        #check if element size is set to points and change meshmode if necessary
+        for pt in self.points:
+            if pt.esize != None:
+                if meshmode=='fineness': print('meshmode is changed to esize, because elementsize was defined on points!')
+                meshmode = 'esize'
+                
+        
+        #if meshmode esize is chosen: ediv's on lines and arcs are transformed to element sizes on start and end point
+        if meshmode == 'esize':
+            for line in self.lines:
+                if line.ediv != None:
+                    line.pt(0).set_esize(line.length()/line.ediv)
+                    line.pt(1).set_esize(line.length()/line.ediv)
+            
+        
         if mesher == 'gmsh':
-            self.__mesh_gmsh(fineness)
+            self.__mesh_gmsh(size, meshmode)
         elif mesher == 'cgx':
-            self.__mesh_cgx(fineness)
+            self.__mesh_cgx(size)
 
-    def __mesh_gmsh(self, fineness):
+    def __mesh_gmsh(self, size, meshmode):
         """Meshes all parts using the Gmsh mesher.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
-
-                Low numbers are very fine, higher numbers are coarser.
+            
+            size (float): 
+            
+                - if meshmode == 'fineness' (default):
+                    - mesh size is adapted to geometry size
+                    - set size = 0.0001 - 1.0, to define how fine the mesh is.
+                    - Low numbers are very fine, higher numbers are coarser.
+                
+                - if meshmode == 'esize':
+                    - element size is kept constant
+                    - choose it depending on geometry size
+                
+            meshmode (str):
+                
+                - 'fineness': adapt mesh size to geometry
+                - 'esize': keep explicitly defined element size
+            
+            
         """
         geo = []
         ids = {}
@@ -1525,6 +1589,15 @@ def __mesh_gmsh(self, fineness):
         # write all points
         for pt in self.points:
             txtline = 'Point(%i) = {%f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0)
+            
+            if meshmode == 'esize':
+                #add element size to points
+                if pt.esize == None:
+                    txtline = txtline.replace('}', ', %f}' % (size if self.eshape=='tri' else size*2.))
+                    #txtline = txtline.replace('}', ', %f}' % (size))
+                else:
+                    txtline = txtline.replace('}', ', %f}' % (pt.esize if self.eshape=='tri' else pt.esize*2.))
+                    #txtline = txtline.replace('}', ', %f}' % (pt.esize))
             geo.append(txtline)
 
         # start storing an index number
@@ -1545,7 +1618,7 @@ def __mesh_gmsh(self, fineness):
             geo.append(txtline)
 
             # set division if we have it
-            if line.ediv != None:
+            if line.ediv != None and meshmode=='fineness':
                 ndiv = line.ediv+1
                 esize = line.length()/line.ediv
                 if self.eshape == 'quad':
@@ -1613,8 +1686,8 @@ def __mesh_gmsh(self, fineness):
             geo.append(txtline)
 
         # set the meshing options
-        geo.append('Mesh.CharacteristicLengthFactor = '
-                   +str(fineness)+'; //mesh fineness')
+        if meshmode == 'fineness':
+            geo.append('Mesh.CharacteristicLengthFactor = '+str(size)+'; //mesh fineness')
         geo.append('Mesh.RecombinationAlgorithm = 1; //blossom')
 
         if self.eshape == 'quad':
@@ -1663,13 +1736,26 @@ def __mesh_gmsh(self, fineness):
         # read in the calculix mesh
         self.__read_inp(self.fname+'.inp')
 
-    def __mesh_cgx(self, fineness):
+    def __mesh_cgx(self, size, meshmode):
         """Meshes all parts using the Calculix cgx mesher.
 
         Args:
-            fineness (float): 0.0001 - 1.0, how fine the mesh is.
-
-                Low numbers are very fine, higher numbers are coarser.
+            
+            size (float): 
+            
+                - if meshmode == 'fineness' (default):
+                    - mesh size is adapted to geometry size
+                    - set size = 0.0001 - 1.0, to define how fine the mesh is.
+                    - Low numbers are very fine, higher numbers are coarser.
+                
+                - if meshmode == 'esize':  NOT TESTED WITH CGX
+                    - element size is kept constant
+                    - choose it depending on geometry size
+                
+            meshmode (str):
+                
+                - 'fineness': adapt mesh size to geometry
+                - 'esize': keep explicitly defined element size  NOT TESTED WITH CGX
         """
         fbd = []
         comps = []
@@ -1693,8 +1779,8 @@ def __mesh_cgx(self, fineness):
         cgx_elements['quad2plstrain'] = 'qu8e'
         cgx_elements['quad1plstrain'] = 'qu4e'
 
-        num = 1.0/fineness
-        emult = int(round(num)) # this converts fineness to mesh multiplier
+        num = 1.0/size
+        emult = int(round(num)) # this converts size to mesh multiplier
 
         # write all points
         for point in self.points:
diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py
index 51e6ed5..749d6f5 100644
--- a/pycalculix/geometry.py
+++ b/pycalculix/geometry.py
@@ -67,6 +67,7 @@ def __init__(self, x, y, z=0):
         self.lines = []
         self.arc_center = False
         self.on_part = True
+        self.esize = None
         base_classes.Idobj.__init__(self)
 
     def get_name(self):
@@ -258,6 +259,10 @@ def __str__(self):
         """Returns string listing object type, id number, and coordinates"""
         val = 'Point %s, (x, y)=(%.3f, %.3f)' % (self.get_name(), self.x, self.y)
         return val
+    
+    def set_esize(self, size):
+        self.esize = size
+        
 
 
 class Line(base_classes.Idobj):
@@ -352,6 +357,15 @@ def set_ediv(self, ediv):
             ediv (int): number of required elements on this line
         """
         self.ediv = ediv
+        
+    def set_esize(self, esize):
+        """Sets the size of mesh elements on the line when meshing.
+
+        Args:
+            esize (float): size of mesh elements on this line
+        """
+        for i in range(2):
+            self.pt(i).set_esize(esize)
 
     def signed_copy(self, sign):
         """Returns a SignLine copy of this Line with the passed sign."""
@@ -675,7 +689,11 @@ def set_pt(self, index, point):
     def set_ediv(self, ediv):
         """Applies the element divisions onto the parent line."""
         self.line.set_ediv(ediv)
-
+        
+    def set_esize(self, esize):
+        """Applies the element size onto the parent line."""
+        self.line.set_esize(esize)
+        
     def set_lineloop(self, lineloop):
         """Sets the parent LineLoop"""
         # this is needed to cascade set ediv up to FEA model and down onto
@@ -869,7 +887,7 @@ def save_to_points(self):
         """This stores this line in the line's child points."""
         for point in self.allpoints:
             point.save_line(self)
-
+        
     def set_ediv(self, ediv):
         """Sets the number of element divisions on the arc when meshing.
 
@@ -877,7 +895,17 @@ def set_ediv(self, ediv):
             ediv (int): number of required elements on this arc
         """
         self.ediv = ediv
+        
+    def set_esize(self, esize):
+        """Sets the size of mesh elements on the arc when meshing.
 
+        Args:
+            esize (float): size of mesh elements on this arc
+        """
+        
+        for i in range(2):
+            self.pt(i).set_esize(esize)
+        
     def signed_copy(self, sign):
         """Returns a SignArc instance of this Arc with the passed sign."""
         return SignArc(self, sign)
@@ -1256,6 +1284,10 @@ def set_pt(self, ind, point):
     def set_ediv(self, ediv):
         """Apply the element divisions onto the parent line."""
         self.line.set_ediv(ediv)
+        
+    def set_esize(self, esize):
+        """Applies the element size onto the parent line."""
+        self.line.set_esize(esize)
 
     def set_lineloop(self, lineloop):
         """Sets the parent LineLoop"""
@@ -1657,7 +1689,7 @@ class Area(base_classes.Idobj):
         signlines (list): list of signed lines or arcs that define the area
         lines (list): a list of all the lines that make the area, includes hole
             lines
-        point (list): a list of all points making the area, excludes arc centers
+        points (list): a list of all points making the area, excludes arc centers
         allpoints (list): a list of all points, includes arc centers
         holepoints (list): a list of hole points, excludes arc centers
         matl (Matl): the material fo the area
@@ -1983,6 +2015,16 @@ def set_etype(self, etype):
         """
         self.etype = etype
         self.set_child_ccxtypes()
+        
+    def set_esize(self, esize):
+        """Sets the size of mesh elements on the area when meshing.
+
+        Args:
+            esize (float): size of mesh elements on this area
+        """
+        #size if set if no mesh size on points is peresent
+        for p in self.points:
+            p.set_esize(esize)
 
     def update(self, line_list):
         """Updates the area's external lines to the passed list.

From 96c06353bcc1d51d49fda0bcd0e82ffadee8b095 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Sun, 5 Jul 2015 18:07:41 +0200
Subject: [PATCH 6/9] add nonlinear mechanical material add material
 attributes: - mechtype = 'linear' or 'nonlinear' - exponent: exponent of
 Ramberg-Osgood stress-strain equation - yield_stress: yield stress of
 material - yield_offset: yield offset of Ramberg-Osgood stress-strain
 equation

Default setting is 'linear' material and even if the nonlinear variables are set, they are not used until mechtpye='nonlinear' is explicitly called.
---
 pycalculix/material.py | 63 +++++++++++++++++++++++++++++-------------
 pycalculix/problem.py  | 10 +++----
 2 files changed, 49 insertions(+), 24 deletions(-)

diff --git a/pycalculix/material.py b/pycalculix/material.py
index 7542117..b39aad9 100644
--- a/pycalculix/material.py
+++ b/pycalculix/material.py
@@ -10,16 +10,22 @@ class Material(base_classes.Idobj):
         name (str): a unique matl name
 
     Attributes:
-        id (int): material id number
-        density (float): density in mass/volume units
-        pratio (float): poisson's ratio (unitless)
-        youngs (float): young's modulus in Force/area = stress units
-        conductivity (float): thermal conductivity, Power / (distance-temp)
-        spec_heat (float): specific heat (energy/(mass-temp)
-        thermal_exp (dict): a dict storing temperature dependent thermal
-            expansion properties
+        - id (int): material id number
+        - mechtype (str): defines mechanical material type: linear or non-linear material
+        - density (float): density in mass/volume units
+        - pratio (float): poisson's ratio (unitless)
+        - youngs (float): young's modulus in Force/area = stress units
+        - mechtpye (str): mechanical material type
+            options: 'linear' or 'nonlinear'
+        - exponent (float): exponent of Ramberg-Osgood stress-strain equation
+        - yield_stress (float): yield stress of material
+        - yield_offset (float): yield offset of Ramberg-Osgood stress-strain equation
+        - conductivity (float): thermal conductivity, Power / (distance-temp)
+        - spec_heat (float): specific heat (energy/(mass-temp)
+        - thermal_exp (dict): a dict storing temperature dependent thermal
+            -- expansion properties
 
-            Thermal expansion is in strain per temperature
+            -- Thermal expansion is in strain per temperature
 
             - dict['data'] = zip of (temp, thermal_expansion)
             - dict['tzero'] = the temperature zero point
@@ -29,33 +35,48 @@ def __init__(self, name):
         self.name = name
         base_classes.Idobj.__init__(self)
         # mechanical
+        self.mechtype = None
         self.density = None
         self.youngs = None
         self.pratio = None
+        self.exponent = None
+        self.yield_stress = None
+        self.yield_offset = None
         # thermal
         self.conductivity = None
         self.spec_heat = None
         # thermal growth
         self.thermal_exp = {}
 
-    def set_mech_props(self, density, youngs, pratio):
+    def set_mech_props(self, density, youngs, pratio, mechtype='linear', exponent=0., yield_stress=0., yield_offset=0.):
         """Sets the mechanical properties: density, youngs, poisson_ratio.
 
         Args:
-          density (float): density in mass/volume units
-          pratio (float): poisson's ratio (unitless)
-          youngs (float): young's modulus in Force/area = stress units
+          - density (float): density in mass/volume units
+          - pratio (float): poisson's ratio (unitless)
+          - youngs (float): young's modulus in Force/area = stress units
+          
+        Kargs:
+          - mechtpye (str): mechanical material type
+            options: 'linear' or 'nonlinear'
+          - exponent (float): exponent of Ramberg-Osgood stress-strain equation
+          - yield_stress (float): yield stress of material
+          - yield_offset (float): yield offset of Ramberg-Osgood stress-strain equation
         """
         self.density = density
         self.youngs = youngs
         self.pratio = pratio
+        self.mechtype = mechtype
+        self.exponent = exponent
+        self.yield_stress = yield_stress
+        self.yield_offset = yield_offset
 
     def set_therm_props(self, conductivity, spec_heat):
         """Sets the thermal properties: conductivity, specifc_heat.
 
         Args:
-          conductivity (float): Power / (distance-temp)
-          spec_heat (float): specific heat (energy/(mass-temp)
+          - conductivity (float): Power / (distance-temp)
+          - spec_heat (float): specific heat (energy/(mass-temp)
         """
         self.conductivity = conductivity
         self.spec_heat = spec_heat
@@ -64,10 +85,10 @@ def set_therm_expan(self, alphas, temps=None, tzero=None):
         """Sets the thermal expansion of the material.
 
         Args:
-          alphas (float or list): list of thermal expansion alphas
+          - alphas (float or list): list of thermal expansion alphas
                          length must be the same as the passed temps
-          temps (list): list of temperatures
-          tzero (float): temperature zero point
+          - temps (list): list of temperatures
+          - tzero (float): temperature zero point
         """
         self.thermal_exp = {}
         isnum = (isinstance(alphas, float) or isinstance(alphas, int))
@@ -82,9 +103,13 @@ def ccx(self):
         """Returns a list of text strings for ccx defining the material."""
         res = []
         res.append('*MATERIAL,NAME='+self.name)
-        if self.youngs != None:
+        
+        if self.mechtype == 'linear':
             res.append('*ELASTIC')
             res.append(str(self.youngs)+','+str(self.pratio))
+        elif self.mechtype == 'nonlinear':
+            res.append('*DEFORMATION PLASTICITY')
+            res.append(','.join([str(self.youngs),str(self.pratio),str(self.yield_stress),str(self.exponent),str(self.yield_offset)]))
         if self.density != None:
             res.append('*DENSITY')
             res.append(str(self.density))
diff --git a/pycalculix/problem.py b/pycalculix/problem.py
index 6f743df..2c1abaf 100644
--- a/pycalculix/problem.py
+++ b/pycalculix/problem.py
@@ -13,12 +13,12 @@ class Problem(base_classes.Idobj):
     """Makes a problem which can be analyzed with Calculix ccx.
 
     Args:
-        feamodel (FeaModel): the parent FeaModel
-        problem_type (str): model type, options:
-
-            - 'struct': structural
-        fname (str): file prefix for the problem .inp and results files
+        - feamodel (FeaModel): the parent FeaModel
+        - problem_type (str): model type, options:
+            -- 'struct': structural
+        - fname (str): file prefix for the problem .inp and results files
             If value is '' it will default to the project name of the FeaModel
+            
     Attributes:
         fea (FeaModel): parent FeaModel
         __ptype (str): problem type, options:

From 44b1504fe942d8982802dd251a18c78ac16f0c90 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Sun, 5 Jul 2015 18:21:55 +0200
Subject: [PATCH 7/9] add gradient functionality

Added function plot_gradient creates plot where field values are projected onto a given line. The line is defined by start and end point. The line is separated into subpoints (number can be chosen by user) and the field value, which is next to the suboints is used. An interpolation, where the distance of the nodes to the subpoints is used to calculate a weighted value didn't worked out so well, therefore this simple approach.
A poly line can be fitted to the resulting gradient and the equation is shown in the plot legend.

Added function get_relative_gradient calculates the relative (e.g. stress) gradient via a fitted poly line. Parts of the above function are used.
---
 pycalculix/results_file.py | 177 +++++++++++++++++++++++++++++++++++++
 1 file changed, 177 insertions(+)

diff --git a/pycalculix/results_file.py b/pycalculix/results_file.py
index d7a62e2..06db0e5 100644
--- a/pycalculix/results_file.py
+++ b/pycalculix/results_file.py
@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 import matplotlib.colors as colors
 import matplotlib.cm as cmx
+import numpy as np
 
 from . import base_classes # needed for RESFIELDS
 from . import mesh
@@ -424,6 +425,182 @@ def set_time(self, time):
         else:
             print('Time %f is not in the loaded times. Valid times are:')
             print(self.steps)
+    
+    
+    def plot_gradient(self, start_point, end_point, field, fname='', display=True, title='', max_val=None, min_val=None, curve_fitting=True, n_poly=3, n_subpoints=500, legend=True):
+        """Create diagram with data projected onto line on the undeformed geometry.
+
+        Args:
+            start_point [(float), (float)]: starting point of line. [x, y]
+            end_point [(float), (float)]: end point of line. Example: [x, y]
+            field (str): results item to plot, examples: 'ux', 'ey', 'Seqv'
+            
+        Kargs:
+            fname (str): prefix of png file name, if writing an image
+            display (bool): True = interactively show the plot
+            title (str): third line in the plot title
+            max_val (float or None): max value in the y-axis
+                - None: max from selected data used
+                - float: use the passed float
+            min_val (float or None): min value in the y-axis
+                - None: min from selected data used
+                - float: use the passed float
+            curve_fitting (bool): True = a curve is fitted to the gradient
+            n_poly (int): numbers of polygons for fitting
+            n_subpoints (int): numbers of points the line is subdivided into
+            legend (bool): True = legend with fitted equation is shown
+        """
+        
+        # store the selected nodes and elements
+        sel = {}
+        sel['nodes'] = self.__problem.fea.view.nodes
+
+        # sort nodes low to high so index is correct
+        # we have index to id below so showing subsets works
+        sel['nodes'] = list(sel['nodes'])
+        sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+        # store results at nodes
+        node_position = np.zeros((len(sel['nodes']),2))
+        field_values = np.zeros(len(sel['nodes']))
+        
+        for idx, node in enumerate(sel['nodes']):
+            
+            node_position[idx] = [node.x, node.y]
+            field_values[idx] = self.__results[self.__time]['node'][node.id][field]
+            
+        
+        #create subpoints on line
+        subpoints = np.zeros((n_subpoints, 3))  #[x, y, line position]
+        
+        subpoints[:,0] = np.linspace(start_point[0], end_point[0], n_subpoints)
+        subpoints[:,1] = np.linspace(start_point[1], end_point[1], n_subpoints)
+        subpoints[:,2] = np.arange(n_subpoints) / n_subpoints * np.sqrt(np.sum( (np.array(start_point) - np.array(end_point))**2))
+        
+        #calculate weighted field value for every subpoint
+        wfield = np.zeros(n_subpoints)
+        
+        for idx in range(n_subpoints):
+            
+            #calculate inverse of distance from nodes to subpoints
+            dist = np.sqrt(np.sum((node_position-subpoints[idx,0:2])**2,axis=1))
+            
+            #calculte weighted field value
+            #dist[dist < 1E-10] = 1E-10
+            #inv_dist = 1. / dist**3
+            #wfield[idx] = np.average(field_values, weights=inv_dist)
+            
+            #use nearest value
+            wfield[idx] = field_values[min(range(len(dist)),key=dist.__getitem__)]
+            
+            
+        #plot diagram
+        
+        fig = plt.figure(figsize=(10,6))
+        ax_ = fig.add_subplot(111)
+        
+        plt.plot(subpoints[:,2], wfield, '-r', linewidth=2.5, label=field)
+        
+        if curve_fitting==True:
+            #execute curve fitting if needed
+            poly = np.polyfit(subpoints[:,2], wfield, n_poly)
+            
+            #string for equation of fitted function
+            funcstring = [str(np.round(poly[i]))+u'*x^'+str(np.arange(n_poly,0,-1)[i]) for i in range(n_poly)]
+            funcstring.append(str(np.round(poly[-1])))
+            funcstring = '+'.join(funcstring)
+            
+            func = np.poly1d(poly)
+            
+            plt.plot(subpoints[:,2], func(subpoints[:,2]), '--k', linewidth=1.5, label=funcstring)
+        
+        
+        # set units
+        alist = self.__problem.fea.get_units(field, 'dist', 'time')
+        [f_unit, d_unit, t_unit] = alist
+
+        # set plot axes
+        plot_title = ('Gradient %s%s\nTime=%f%s' %(field, f_unit, self.__time, t_unit))
+        if title != '':
+            plot_title += '\n%s' % title
+        plt.title(plot_title)
+        plt.xlabel('path position'+d_unit)
+        plt.ylabel(field + ' ' +f_unit)
+        
+        #show legend if needed
+        if legend == True:
+            plt.legend()
+            
+        #set limits on y-axis
+        if min_val!=None:
+            plt.gca().set_ylim(bottom=min_val)
+        if max_val!=None:
+            plt.gca().set_ylim(top=max_val)
+            
+        plt.grid()
+        base_classes.plot_finish(plt, fname, display)
+    
+    def get_relative_gradient(self, start_point, end_point, field, n_poly=3, n_subpoints=500):
+        """Calculte relative stress gradient (gradient/start_value)
+
+        Args:
+            start_point [(float), (float)]: starting point of line. [x, y]
+            end_point [(float), (float)]: end point of line. Example: [x, y]
+            field (str): results item to plot, examples: 'ux', 'ey', 'Seqv'
+            
+        Kargs:
+            n_poly (int): numbers of polygons for fitting, min=2
+            n_subpoints (int): numbers of points the line is subdivided into
+        """
+        
+        # store the selected nodes and elements
+        sel = {}
+        sel['nodes'] = self.__problem.fea.view.nodes
+
+        # sort nodes low to high so index is correct
+        # we have index to id below so showing subsets works
+        sel['nodes'] = list(sel['nodes'])
+        sel['nodes'] = sorted(sel['nodes'], key=lambda k: k.id)
+
+        # store results at nodes
+        node_position = np.zeros((len(sel['nodes']),2))
+        field_values = np.zeros(len(sel['nodes']))
+        
+        for idx, node in enumerate(sel['nodes']):
+            
+            node_position[idx] = [node.x, node.y]
+            field_values[idx] = self.__results[self.__time]['node'][node.id][field]
+            
+        
+        #create subpoints on line
+        subpoints = np.zeros((n_subpoints, 3))  #[x, y, line position]
+        
+        subpoints[:,0] = np.linspace(start_point[0], end_point[0], n_subpoints)
+        subpoints[:,1] = np.linspace(start_point[1], end_point[1], n_subpoints)
+        subpoints[:,2] = np.arange(n_subpoints) / n_subpoints * np.sqrt(np.sum( (np.array(start_point) - np.array(end_point))**2))
+        
+        #calculate weighted field value for every subpoint
+        wfield = np.zeros(n_subpoints)
+        
+        for idx in range(n_subpoints):
+            
+            #calculate inverse of distance from nodes to subpoints
+            dist = np.sqrt(np.sum((node_position-subpoints[idx,0:2])**2,axis=1))
+            
+            #use nearest value
+            wfield[idx] = field_values[min(range(len(dist)),key=dist.__getitem__)]
+            
+            
+        #curve fitting
+        poly = np.polyfit(subpoints[:,2], wfield, n_poly)
+        
+        rel_grad = abs(poly[-2])/abs(poly[-1])
+        
+        return rel_grad
+        
+        
+        
+        
 
     @staticmethod
     def __utot(vals):

From 6a77a489770ffe1ff905af53b2d95ce298fa2385 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Sun, 5 Jul 2015 18:22:24 +0200
Subject: [PATCH 8/9] Revert "Fixed meshing problems with gmsh"

This reverts commit bd1bb506f07d543df0c7689bab342db7c1914501.
---
 pycalculix/geometry.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py
index eb53522..8bc20cb 100644
--- a/pycalculix/geometry.py
+++ b/pycalculix/geometry.py
@@ -356,7 +356,7 @@ def set_ediv(self, ediv):
         Args:
             ediv (int): number of required elements on this line
         """
-        
+        print('Funktion ediv aufgerufen')
         self.ediv = ediv
         self.pt(0).set_esize(self.length()/ediv)
         self.pt(1).set_esize(self.length()/ediv)
@@ -893,20 +893,20 @@ def save_to_points(self):
             point.save_line(self)
         
     def set_ediv(self, ediv):
-        """Sets the number of element divisions on the arc when meshing.
+        """Sets the number of element divisions on the line when meshing.
 
         Args:
-            ediv (int): number of required elements on this arc
+            ediv (int): number of required elements on this line
         """
         self.ediv = ediv
         self.pt(0).set_esize(self.length()/ediv)
         self.pt(1).set_esize(self.length()/ediv)
         
     def set_esize(self, esize):
-        """Sets the size of mesh elements on the arc when meshing.
+        """Sets the size of mesh elements on the line when meshing.
 
         Args:
-            esize (float): size of mesh elements on this arc
+            esize (float): size of mesh elements on this line
         """
         self.pt(0).set_esize(esize)
         self.pt(1).set_esize(esize)

From 4be77f45bbef77315a2338dd9e8b3d6b26b3c2e8 Mon Sep 17 00:00:00 2001
From: KyPy <kai.sandmann@web.de>
Date: Sun, 5 Jul 2015 18:28:28 +0200
Subject: [PATCH 9/9] Undo prior changes synchronizing branch with main branch
 of spacether

---
 pycalculix/connectors.py |  14 ++--
 pycalculix/feamodel.py   | 155 +++++++++++++++++++++++++++++----------
 pycalculix/geometry.py   |  41 +----------
 pycalculix/mesh.py       |   1 +
 4 files changed, 128 insertions(+), 83 deletions(-)

diff --git a/pycalculix/connectors.py b/pycalculix/connectors.py
index b33c0f6..e4f077f 100644
--- a/pycalculix/connectors.py
+++ b/pycalculix/connectors.py
@@ -10,7 +10,7 @@ class SurfaceInteraction(base_classes.Idobj):
 
             * 'EXPONENTIAL'
             * 'LINEAR'
-        
+
         *args: following arguments
 
             * int_type = 'EXPONENTIAL'
@@ -20,7 +20,7 @@ class SurfaceInteraction(base_classes.Idobj):
                 * k must be passed
 
     """
-    
+
     def __init__(self, int_type, *args):
         self.int_type = int_type
         self.c0 = None
@@ -32,9 +32,10 @@ def __init__(self, int_type, *args):
         elif self.int_type == 'LINEAR':
             self.k = args[0]
         base_classes.Idobj.__init__(self)
-    
+
     @property
     def name(self):
+        """SurfaceInteraction name."""
         return 'SI%i' % self.id
 
     def ccx(self):
@@ -50,7 +51,7 @@ def ccx(self):
 
 class Contact(base_classes.Idobj):
     """Makes a contact which will be between lines which have faces on them.
-    
+
     Args:
         master_comp (Component): component of master element faces
         slave_comp (Component): component of slave element faces
@@ -58,7 +59,7 @@ class Contact(base_classes.Idobj):
         surf_to_surf (bool): if True surface to surface is used, if False node
             to surface is used
     """
-    
+
     def __init__(self, master_comp, slave_comp, surf_int, surf_to_surf=True):
         self.master_comp = master_comp
         self.slave_comp = slave_comp
@@ -68,6 +69,7 @@ def __init__(self, master_comp, slave_comp, surf_int, surf_to_surf=True):
 
     @property
     def name(self):
+        """Contact name."""
         return 'CONT%i' % self.id
 
     def ccx(self):
@@ -83,4 +85,4 @@ def ccx(self):
         res.append(line)
         line = "%s,%s" % (s_name, m_name)
         res.append(line)
-        return res
\ No newline at end of file
+        return res
diff --git a/pycalculix/feamodel.py b/pycalculix/feamodel.py
index 53cee80..85adf66 100644
--- a/pycalculix/feamodel.py
+++ b/pycalculix/feamodel.py
@@ -67,6 +67,8 @@ class FeaModel(object):
             | Time = 0.0 stores constant loads, such as:
             |    material, thickness
 
+        contacts (Itemlist): list of contacts
+        surfints (Itemlist): list of surface interactions
         problems (Itemlist): list of problems
         nodes (Meshlist): list of all mesh nodes
         eshape (str): element shape
@@ -158,22 +160,6 @@ def set_ediv(self, items, ediv):
         items = self.get_items(items)
         for line in items:
             line.set_ediv(ediv)
-            
-    def set_esize(self, items, esize):
-        """Sets the element size on the passed line.
-
-        Args:
-            items (str or SignLine or SignArc or list): lines or points to set esize on
-
-                - str: 'L0'
-                - list of str ['L0', 'L1', 'P3']
-                - list of SignLine or SignArc part.bottom or part.hole[-1]
-
-            esize (float): size of the mesh elements on the line
-        """
-        items = self.get_items(items)
-        for item in items:
-            item.set_esize(esize)
 
     def set_units(self, dist_unit='m', cfswitch=False):
         """Sets the units that will be displayed when plotting.
@@ -367,6 +353,54 @@ def print_summary(self):
             print(' %s: %i' % (name, len(items)))
         print(spacer)
 
+    def plot_nodes(self, fname='', display=True, title='Nodes', nnum=False):
+        """Plots the selected nodes.
+
+        Args:
+            fname (str): png image file prefix, if given, image will be saved
+            display (bool): if True, interactive plot will be shown, if False
+                plot will not be shown
+            title (str): the plot's title
+            nnum (bool): if True node numbers are plotted
+        """
+        nodes = self.view.nodes
+        if len(nodes) > 0:
+            # plotting elements
+            fig = plt.figure()
+            axis = fig.add_subplot(111)
+
+            # plot nodes, this is quicker than individual plotting
+            axs = [node.y for node in nodes]
+            rads = [node.x for node in nodes]
+            axis.scatter(axs, rads, s=7, color='black')
+            if nnum:
+                for node in nodes:
+                    node.label(axis)
+
+            # set units
+            [d_unit] = self.get_units('dist')
+            plt.title(title)
+            plt.xlabel('axial, y'+d_unit)
+            plt.ylabel('radial, x'+d_unit)
+            plt.axis('scaled')
+
+            # extract max and min for plot window
+            radials = [n.x for n in nodes]
+            axials = [n.y for n in nodes]
+
+            # finish pot
+            base_classes.plot_set_bounds(plt, axials, radials)
+            base_classes.plot_finish(plt, fname, display)
+
+        else:
+            # no elements exist or no elemnts are selected
+            res = ''
+            if len(self.nodes) == 0:
+                res = 'No nodes exist! Try meshing your parts!'
+            else:
+                res = 'No nodes are selected! Select some!'
+            print(res)
+
     def plot_elements(self, fname='', display=True, title='Elements',
                       enum=False, nshow=False, nnum=False):
         """Plots the selected elements.
@@ -559,7 +593,7 @@ def plot_pressures(self, fname='', display=True):
             cbar = plt.colorbar(scalarmap, orientation='vertical',
                                 ticks=tick_list)
             if cbar_val != None:
-                cbar.ax.set_yticklabels(['',str(cbar_val),''])
+                cbar.ax.set_yticklabels(['', str(cbar_val), ''])
             base_classes.plot_finish(plt, fname, display)
 
         else:
@@ -688,7 +722,7 @@ def plot_constraints(self, fname='', display=True):
             cbar = plt.colorbar(scalarmap, orientation='vertical',
                                 ticks=tick_list)
             if cbar_val != None:
-                cbar.ax.set_yticklabels(['',str(cbar_val),''])
+                cbar.ax.set_yticklabels(['', str(cbar_val), ''])
             base_classes.plot_finish(plt, fname, display)
 
         else:
@@ -918,6 +952,24 @@ def __get_make_comp(self, comp):
             comp = self.components[ind]
         return comp
 
+    def __get_make_surfint(self, surfint):
+        """Stores surfac interaction if it doesn't exist, returns it if it does.
+
+        Args:
+            surfint (connectors.SurfaceInteraction): item to get or make
+        """
+        items = [item for item in self.surfints if item.int_type == surfint.int_type]
+        if surfint.int_type == 'LINEAR':
+            for item in items:
+                if item.k == surfint.k:
+                    return item
+        elif surfint.type == 'EXPONENTIAL':
+            for item in items:
+                if item.c0 == surfint.c0 and item.p0 == surfint.p0:
+                    return item
+        surfint = self.surfints.append(surfint)
+        return surfint
+
     def register(self, item):
         """Adds an item to the feamodel.
 
@@ -1203,12 +1255,18 @@ def set_constr(self, ltype, items, axis, val=0.0):
         self.__add_load(load, self.time)
         return load
 
-    def set_contact_linear(self, master_lines, slave_lines, kval):
+    def set_contact_linear(self, master_lines, slave_lines, kval, many_si=False):
         """Sets contact between master and slave lines.
-        
+
+        Slave lines are on the more compliant or more curved object.
+
         Args:
             master_lines (list): list of SignLine or SignArc
             slave_lines (list): list of SignLine or SignArc
+            kval (float): stiffness, 5 to 50 times the youngs modulus
+                of the touching matl
+            many_si (bool): True, make unique surface interaction for every contact
+                False, use existing surface interaction if we can
         """
 
         master_items = self.get_items(master_lines)
@@ -1222,10 +1280,13 @@ def set_contact_linear(self, master_lines, slave_lines, kval):
         ctype = 'faces'
         slave_comp = components.Component(slave_items, ctype, slave_cname)
         slave_comp = self.__get_make_comp(slave_comp)
-        
+
         surf_int = connectors.SurfaceInteraction('LINEAR', kval)
-        self.surfints.append(surf_int)
-        
+        if many_si:
+            surf_int = self.surfints.append(surf_int)
+        else:
+            surf_int = self.__get_make_surfint(surf_int)
+
         cont = connectors.Contact(master_comp, slave_comp, surf_int, True)
         self.contacts.append(cont)
 
@@ -1363,7 +1424,7 @@ def __read_inp(self, fname):
                     L = line.split(',')
                     L = [int(a.strip()) for a in L]
                     enum = L[0]
-                    nlist = [Dict_NodeIDs[a] for a in L[1:]]
+                    nlist = [N.idget(a) for a in L[1:]]
                     e = mesh.Element(enum, etype, nlist)
                     faces = e.faces
                     E.append(e)
@@ -1416,7 +1477,7 @@ def __read_inp(self, fname):
                 sets[set_type][set_name] = []
                 mode = 'set'
         f.close()
-        
+
         # loop through sets and remove empty sets
         # store sets to delete
         todel = []
@@ -1441,7 +1502,7 @@ def __read_inp(self, fname):
         self.elements = E
         self.faces = F
 
-        # remove arc center ndoes from imported node set
+        # remove arc center nodes from imported node set
         # those nodes have no elements under them
         torem = []
         for node in N:
@@ -1508,11 +1569,11 @@ def __read_inp(self, fname):
         print('Nodes: %i' % len(N))
         print('Done reading Calculix/Abaqus .inp file')
 
-    def mesh(self, element_size=1.0, mesher='gmsh'):
+    def mesh(self, fineness=1.0, mesher='gmsh'):
         """Meshes all parts.
 
         Args:
-            element_size (float): 0.0001 - 1.0, how fine the mesh is.
+            fineness (float): 0.0001 - 1.0, how fine the mesh is.
 
                 - Low numbers are very fine, higher numbers are coarser.
             mesher (str): the mesher to use
@@ -1521,15 +1582,15 @@ def mesh(self, element_size=1.0, mesher='gmsh'):
                 - 'cgx': mesh with Calculix cgx, it doesn't allow holes
         """
         if mesher == 'gmsh':
-            self.__mesh_gmsh(element_size)
+            self.__mesh_gmsh(fineness)
         elif mesher == 'cgx':
-            self.__mesh_cgx(element_size)
+            self.__mesh_cgx(fineness)
 
-    def __mesh_gmsh(self, element_size):
+    def __mesh_gmsh(self, fineness):
         """Meshes all parts using the Gmsh mesher.
 
         Args:
-            element_size (float): 0.0001 - 1.0, how fine the mesh is.
+            fineness (float): 0.0001 - 1.0, how fine the mesh is.
 
                 Low numbers are very fine, higher numbers are coarser.
         """
@@ -1540,10 +1601,7 @@ def __mesh_gmsh(self, element_size):
 
         # write all points
         for pt in self.points:
-            if pt.esize == None:
-                txtline = 'Point(%i) = {%f, %f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0, element_size if self.eshape=='tri' else element_size*2.)
-            else:
-                txtline = 'Point(%i) = {%f, %f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0, pt.esize if self.eshape=='tri' else pt.esize*2.)
+            txtline = 'Point(%i) = {%f, %f, %f};' % (pt.id, pt.x, pt.y, 0.0)
             geo.append(txtline)
 
         # start storing an index number
@@ -1563,6 +1621,21 @@ def __mesh_gmsh(self, element_size):
                 txtline = 'Line(%i) = {%i,%i};' % (ind, pt1, pt2)
             geo.append(txtline)
 
+            # set division if we have it
+            if line.ediv != None:
+                ndiv = line.ediv+1
+                esize = line.length()/line.ediv
+                if self.eshape == 'quad':
+                    ndiv = line.ediv/2+1
+                    esize = esize*2
+                    # this is needed because quad recombine
+                    # splits 1 element into 2
+                txtline = 'Transfinite Line{%i} = %i;' % (ind, ndiv)
+                print('LINE ELEMENT SIZE: %f, MAKES %i ELEMENTS'
+                      % (line.length()/line.ediv, line.ediv))
+                geo.append(txtline)
+                geo.append('Characteristic Length {%i,%i} = %f;'
+                           % (pt1, pt2, esize))
             ind += 1
 
         # write all areas
@@ -1617,6 +1690,8 @@ def __mesh_gmsh(self, element_size):
             geo.append(txtline)
 
         # set the meshing options
+        geo.append('Mesh.CharacteristicLengthFactor = '
+                   +str(fineness)+'; //mesh fineness')
         geo.append('Mesh.RecombinationAlgorithm = 1; //blossom')
 
         if self.eshape == 'quad':
@@ -1665,11 +1740,11 @@ def __mesh_gmsh(self, element_size):
         # read in the calculix mesh
         self.__read_inp(self.fname+'.inp')
 
-    def __mesh_cgx(self, element_size):
+    def __mesh_cgx(self, fineness):
         """Meshes all parts using the Calculix cgx mesher.
 
         Args:
-            element_size (float): 0.0001 - 1.0, how fine the mesh is.
+            fineness (float): 0.0001 - 1.0, how fine the mesh is.
 
                 Low numbers are very fine, higher numbers are coarser.
         """
@@ -1695,8 +1770,8 @@ def __mesh_cgx(self, element_size):
         cgx_elements['quad2plstrain'] = 'qu8e'
         cgx_elements['quad1plstrain'] = 'qu4e'
 
-        num = 1.0/element_size
-        emult = int(round(num)) # this converts element_size to mesh multiplier
+        num = 1.0/fineness
+        emult = int(round(num)) # this converts fineness to mesh multiplier
 
         # write all points
         for point in self.points:
diff --git a/pycalculix/geometry.py b/pycalculix/geometry.py
index 8bc20cb..51e6ed5 100644
--- a/pycalculix/geometry.py
+++ b/pycalculix/geometry.py
@@ -67,7 +67,6 @@ def __init__(self, x, y, z=0):
         self.lines = []
         self.arc_center = False
         self.on_part = True
-        self.esize = None
         base_classes.Idobj.__init__(self)
 
     def get_name(self):
@@ -259,10 +258,6 @@ def __str__(self):
         """Returns string listing object type, id number, and coordinates"""
         val = 'Point %s, (x, y)=(%.3f, %.3f)' % (self.get_name(), self.x, self.y)
         return val
-    
-    def set_esize(self, size):
-        self.esize = size
-        
 
 
 class Line(base_classes.Idobj):
@@ -356,20 +351,7 @@ def set_ediv(self, ediv):
         Args:
             ediv (int): number of required elements on this line
         """
-        print('Funktion ediv aufgerufen')
         self.ediv = ediv
-        self.pt(0).set_esize(self.length()/ediv)
-        self.pt(1).set_esize(self.length()/ediv)
-        
-        
-    def set_esize(self, esize):
-        """Sets the size of mesh elements on the line when meshing.
-
-        Args:
-            esize (float): size of mesh elements on this line
-        """
-        self.pt(0).set_esize(esize)
-        self.pt(1).set_esize(esize)
 
     def signed_copy(self, sign):
         """Returns a SignLine copy of this Line with the passed sign."""
@@ -693,11 +675,7 @@ def set_pt(self, index, point):
     def set_ediv(self, ediv):
         """Applies the element divisions onto the parent line."""
         self.line.set_ediv(ediv)
-        
-    def set_esize(self, esize):
-        """Applies the element size onto the parent line."""
-        self.line.set_esize(esize)
-        
+
     def set_lineloop(self, lineloop):
         """Sets the parent LineLoop"""
         # this is needed to cascade set ediv up to FEA model and down onto
@@ -891,26 +869,15 @@ def save_to_points(self):
         """This stores this line in the line's child points."""
         for point in self.allpoints:
             point.save_line(self)
-        
+
     def set_ediv(self, ediv):
-        """Sets the number of element divisions on the line when meshing.
+        """Sets the number of element divisions on the arc when meshing.
 
         Args:
-            ediv (int): number of required elements on this line
+            ediv (int): number of required elements on this arc
         """
         self.ediv = ediv
-        self.pt(0).set_esize(self.length()/ediv)
-        self.pt(1).set_esize(self.length()/ediv)
-        
-    def set_esize(self, esize):
-        """Sets the size of mesh elements on the line when meshing.
 
-        Args:
-            esize (float): size of mesh elements on this line
-        """
-        self.pt(0).set_esize(esize)
-        self.pt(1).set_esize(esize)
-        
     def signed_copy(self, sign):
         """Returns a SignArc instance of this Arc with the passed sign."""
         return SignArc(self, sign)
diff --git a/pycalculix/mesh.py b/pycalculix/mesh.py
index 7296363..6d30f87 100644
--- a/pycalculix/mesh.py
+++ b/pycalculix/mesh.py
@@ -319,6 +319,7 @@ def set_nmid(self, nmid):
         """
         self.nmid = nmid
         nmid.add_face(self)
+        self.nodes = [self.nodes[0], self.nmid, self.nodes[-1]]
 
     def __eq__(self, other):
         """Compare this face to other face. True if node lists equal.