FEM: add module FemMeshTools and move methods in there

src/Mod/Fem/App/CMakeLists.txt

@@ -104,6 +104,7 @@ SET(FemScripts_SRCS
     FemAnalysis.py
     FemCommands.py
     FemBeamSection.py
+    FemMeshTools.py
     FemShellThickness.py
     FemSolverCalculix.py
     FemTools.py

src/Mod/Fem/CMakeLists.txt

@@ -15,6 +15,7 @@ INSTALL(
 
         SelectionObserverFem.py
         TestFem.py
+        FemMeshTools.py
 
         ccxDatReader.py
         ccxFrdReader.py

src/Mod/Fem/FemMeshTools.py

@@ -0,0 +1,631 @@
+# ***************************************************************************
+# *                                                                         *
+# *   Copyright (c) 2016 - Bernd Hahnebach <bernd@bimstatik.org>            *
+# *                                                                         *
+# *   This program is free software; you can redistribute it and/or modify  *
+# *   it under the terms of the GNU Lesser General Public License (LGPL)    *
+# *   as published by the Free Software Foundation; either version 2 of     *
+# *   the License, or (at your option) any later version.                   *
+# *   for detail see the LICENCE text file.                                 *
+# *                                                                         *
+# *   This program is distributed in the hope that it will be useful,       *
+# *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+# *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+# *   GNU Library General Public License for more details.                  *
+# *                                                                         *
+# *   You should have received a copy of the GNU Library General Public     *
+# *   License along with this program; if not, write to the Free Software   *
+# *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  *
+# *   USA                                                                   *
+# *                                                                         *
+# ***************************************************************************
+
+
+__title__ = "Tools for the work with FEM meshes"
+__author__ = "Bernd Hahnebach"
+__url__ = "http://www.freecadweb.org"
+
+
+import FreeCAD
+
+
+def get_femelements_by_references(femmesh, femelement_table, references):
+    '''get the femelements for a list of references
+    '''
+    references_femelements = []
+    for ref in references:
+        ref_femnodes = get_femnodes_by_refshape(femmesh, ref)  # femnodes for the current ref
+        references_femelements += get_femelements_by_femnodes(femelement_table, ref_femnodes)  # femelements for all references
+    return references_femelements
+
+
+def get_femnodes_by_references(femmesh, references):
+    '''get the femnodes for a list of references
+    '''
+    references_femnodes = []
+    for ref in references:
+        references_femnodes += get_femnodes_by_refshape(femmesh, ref)
+
+    # return references_femnodes  # keeps duplicate nodes, keeps node order
+
+    # if nodes are used for nodesets, duplicats should be removed
+    return list(set(references_femnodes))  # removes duplicate nodes, sortes node order
+
+
+def get_femnodes_by_refshape(femmesh, ref):
+    if ref[1]:
+        r = ref[0].Shape.getElement(ref[1])  # Vertex, Edge, Face
+    else:
+        r = ref[0].Shape  # solid
+    # print('  ReferenceShape : ', r.ShapeType, ', ', ref[0].Name, ', ', ref[0].Label, ' --> ', ref[1])
+    if r.ShapeType == 'Vertex':
+        nodes = femmesh.getNodesByVertex(r)
+    elif r.ShapeType == 'Edge':
+        nodes = femmesh.getNodesByEdge(r)
+    elif r.ShapeType == 'Face':
+        nodes = femmesh.getNodesByFace(r)
+    elif r.ShapeType == 'Solid':
+        nodes = femmesh.getNodesBySolid(r)
+    else:
+        print('  No Vertice, Edge, Face or Solid as reference shapes!')
+    return nodes
+
+
+def get_femelement_table(femmesh):
+    """ get_femelement_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }"""
+    femelement_table = {}
+    if is_solid_femmesh(femmesh):
+        for i in femmesh.Volumes:
+            femelement_table[i] = femmesh.getElementNodes(i)
+    elif is_face_femmesh(femmesh):
+        for i in femmesh.Faces:
+            femelement_table[i] = femmesh.getElementNodes(i)
+    elif is_edge_femmesh(femmesh):
+        for i in femmesh.Edges:
+            femelement_table[i] = femmesh.getElementNodes(i)
+    else:
+        FreeCAD.Console.PrintError('Neither solid nor face nor edge femmesh!\n')
+    return femelement_table
+
+
+def get_femelements_by_femnodes(femelement_table, node_list):
+    '''for every femelement of femelement_table
+    if all nodes of the femelement are in node_list,
+    the femelement is added to the list which is returned
+    e: elementlist
+    nodes: nodelist '''
+    e = []  # elementlist
+    for elementID in sorted(femelement_table):
+        nodecount = 0
+        for nodeID in femelement_table[elementID]:
+            if nodeID in node_list:
+                nodecount = nodecount + 1
+        if nodecount == len(femelement_table[elementID]):   # all nodes of the element are in the node_list!
+            e.append(elementID)
+    return e
+
+
+def get_femelement_sets(femmesh, femelement_table, fem_objects):  # fem_objects = FreeCAD FEM document objects
+    # get femelements for reference shapes of each obj.References
+    count_femelements = 0
+    referenced_femelements = []
+    has_remaining_femelements = None
+    for fem_object_i, fem_object in enumerate(fem_objects):
+        obj = fem_object['Object']
+        fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i)  # unique short identifier
+        if obj.References:
+            ref_shape_femelements = []
+            ref_shape_femelements = get_femelements_by_references(femmesh, femelement_table, obj.References)
+            referenced_femelements += ref_shape_femelements
+            count_femelements += len(ref_shape_femelements)
+            fem_object['FEMElements'] = ref_shape_femelements
+        else:
+            has_remaining_femelements = obj.Name
+    # get remaining femelements for the fem_objects
+    if has_remaining_femelements:
+        remaining_femelements = []
+        for elemid in femelement_table:
+            if elemid not in referenced_femelements:
+                remaining_femelements.append(elemid)
+        count_femelements += len(remaining_femelements)
+        for fem_object in fem_objects:
+            obj = fem_object['Object']
+            if obj.Name == has_remaining_femelements:
+                fem_object['FEMElements'] = sorted(remaining_femelements)
+    # check if all worked out well
+    if not femelements_count_ok(femelement_table, count_femelements):
+        FreeCAD.Console.PrintError('Error in get_femelement_sets -- > femelements_count_ok() failed!\n')
+
+
+def get_elset_short_name(obj, i):
+    if hasattr(obj, "Proxy") and obj.Proxy.Type == 'MechanicalMaterial':
+        return 'Mat' + str(i)
+    elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemBeamSection':
+        return 'Beam' + str(i)
+    elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemShellThickness':
+        return 'Shell' + str(i)
+    else:
+        print('Error: ', obj.Name, ' --> ', obj.Proxy.Type)
+
+
+def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj):
+    # force_obj_node_load_table = [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)]
+    force_obj_node_load_table = []
+    sum_ref_edge_length = 0
+    sum_ref_edge_node_length = 0  # for debugging
+    sum_node_load = 0  # for debugging
+    for o, elem in frc_obj.References:
+        elem_o = o.Shape.getElement(elem)
+        sum_ref_edge_length += elem_o.Length
+    if sum_ref_edge_length != 0:
+        force_per_sum_ref_edge_length = frc_obj.Force / sum_ref_edge_length
+    for o, elem in frc_obj.References:
+        elem_o = o.Shape.getElement(elem)
+        ref_edge = elem_o
+
+        # edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
+        edge_table = get_ref_edgenodes_table(femmesh, femelement_table, ref_edge)
+
+        # node_length_table = [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
+        node_length_table = get_ref_edgenodes_lengths(femnodes_mesh, edge_table)
+
+        # node_sum_length_table = { nodeID : Length, ... , nodeID : Length }  LengthSum for each node, one entry for each node
+        node_sum_length_table = get_ref_shape_node_sum_geom_table(node_length_table)
+
+        # node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad }  NodeLoad for each node, one entry for each node
+        node_load_table = {}
+        sum_node_lengths = 0  # for debugging
+        for node in node_sum_length_table:
+            sum_node_lengths += node_sum_length_table[node]  # for debugging
+            node_load_table[node] = node_sum_length_table[node] * force_per_sum_ref_edge_length
+        ratio_refedge_lengths = sum_node_lengths / elem_o.Length
+        if ratio_refedge_lengths < 0.99 or ratio_refedge_lengths > 1.01:
+            FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
+            print('  sum_node_lengths:', sum_node_lengths)
+            print('  refedge_length:  ', elem_o.Length)
+            bad_refedge = elem_o
+        sum_ref_edge_node_length += sum_node_lengths
+
+        elem_info_string = 'node loads on shape: ' + o.Name + ':' + elem
+        force_obj_node_load_table.append((elem_info_string, node_load_table))
+
+    for ref_shape in force_obj_node_load_table:
+        for node in ref_shape[1]:
+            sum_node_load += ref_shape[1][node]  # for debugging
+
+    ratio = sum_node_load / frc_obj.Force
+    if ratio < 0.99 or ratio > 1.01:
+        print('Deviation  sum_node_load to frc_obj.Force is more than 1% :  ', ratio)
+        print('  sum_ref_edge_node_length: ', sum_ref_edge_node_length)
+        print('  sum_ref_edge_length:      ', sum_ref_edge_length)
+        print('  sum_node_load:          ', sum_node_load)
+        print('  frc_obj.Force:          ', frc_obj.Force)
+        print('  the reason could be simply a circle length --> see method get_ref_edge_node_lengths')
+        print('  the reason could also be an problem in retrieving the ref_edge_node_length')
+
+        # try debugging of the last bad refedge
+        print('DEBUGGING')
+        print(bad_refedge)
+
+        print('bad_refedge_nodes')
+        bad_refedge_nodes = femmesh.getNodesByEdge(bad_refedge)
+        print(len(bad_refedge_nodes))
+        print(bad_refedge_nodes)
+        # import FreeCADGui
+        # FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_refedge_nodes
+
+        print('bad_edge_table')
+        # bad_edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
+        bad_edge_table = get_ref_edgenodes_table(femmesh, femelement_table, bad_refedge)
+        print(len(bad_edge_table))
+        bad_edge_table_nodes = []
+        for elem in bad_edge_table:
+            print(elem, ' --> ', bad_edge_table[elem])
+            for node in bad_edge_table[elem]:
+                if node not in bad_edge_table_nodes:
+                    bad_edge_table_nodes.append(node)
+        print('sorted(bad_edge_table_nodes)')
+        print(sorted(bad_edge_table_nodes))   # should be == bad_refedge_nodes
+        # import FreeCADGui
+        # FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_edge_table_nodes
+        # bad_node_length_table = [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
+
+        print('good_edge_table')
+        good_edge_table = delete_duplicate_mesh_elements(bad_edge_table)
+        for elem in good_edge_table:
+            print(elem, ' --> ', bad_edge_table[elem])
+
+        print('bad_node_length_table')
+        bad_node_length_table = get_ref_edgenodes_lengths(femnodes_mesh, bad_edge_table)
+        for n, l in bad_node_length_table:
+            print(n, ' --> ', l)
+
+    return force_obj_node_load_table
+
+
+def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj):
+    # force_obj_node_load_table = [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)]
+    force_obj_node_load_table = []
+    sum_ref_face_area = 0
+    sum_ref_face_node_area = 0  # for debugging
+    sum_node_load = 0  # for debugging
+    for o, elem in frc_obj.References:
+        elem_o = o.Shape.getElement(elem)
+        sum_ref_face_area += elem_o.Area
+    if sum_ref_face_area != 0:
+        force_per_sum_ref_face_area = frc_obj.Force / sum_ref_face_area
+    for o, elem in frc_obj.References:
+        elem_o = o.Shape.getElement(elem)
+        ref_face = elem_o
+
+        # face_table = { meshfaceID : ( nodeID, ... , nodeID ) }
+        face_table = get_ref_facenodes_table(femmesh, femelement_table, ref_face)
+
+        # node_area_table = [ (nodeID, Area), ... , (nodeID, Area) ]  some nodes will have more than one entry
+        node_area_table = get_ref_facenodes_areas(femnodes_mesh, face_table)
+
+        # node_sum_area_table = { nodeID : Area, ... , nodeID : Area }  AreaSum for each node, one entry for each node
+        node_sum_area_table = get_ref_shape_node_sum_geom_table(node_area_table)
+
+        # node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad }  NodeLoad for each node, one entry for each node
+        node_load_table = {}
+        sum_node_areas = 0  # for debugging
+        for node in node_sum_area_table:
+            sum_node_areas += node_sum_area_table[node]  # for debugging
+            node_load_table[node] = node_sum_area_table[node] * force_per_sum_ref_face_area
+        ratio_refface_areas = sum_node_areas / elem_o.Area
+        if ratio_refface_areas < 0.99 or ratio_refface_areas > 1.01:
+            FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
+            print('  sum_node_lengths:', sum_node_areas)
+            print('  refedge_length:  ', elem_o.Area)
+        sum_ref_face_node_area += sum_node_areas
+
+        elem_info_string = 'node loads on shape: ' + o.Name + ':' + elem
+        force_obj_node_load_table.append((elem_info_string, node_load_table))
+
+    for ref_shape in force_obj_node_load_table:
+        for node in ref_shape[1]:
+            sum_node_load += ref_shape[1][node]  # for debugging
+
+    ratio = sum_node_load / frc_obj.Force
+    if ratio < 0.99 or ratio > 1.01:
+        print('Deviation  sum_node_load to frc_obj.Force is more than 1% :  ', ratio)
+        print('  sum_ref_face_node_area: ', sum_ref_face_node_area)
+        print('  sum_ref_face_area:      ', sum_ref_face_area)
+        print('  sum_node_load:          ', sum_node_load)
+        print('  frc_obj.Force:          ', frc_obj.Force)
+        print('  the reason could be simply a circle area --> see method get_ref_face_node_areas')
+        print('  the reason could also be an problem in retrieving the ref_face_node_area')
+
+    return force_obj_node_load_table
+
+
+def get_ref_edgenodes_table(femmesh, femelement_table, refedge):
+    edge_table = {}  # { meshedgeID : ( nodeID, ... , nodeID ) }
+    refedge_nodes = femmesh.getNodesByEdge(refedge)
+    if is_solid_femmesh(femmesh):
+        refedge_fem_volumeelements = []
+        # if at least two nodes of a femvolumeelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
+        for elem in femelement_table:
+            nodecount = 0
+            for node in femelement_table[elem]:
+                if node in refedge_nodes:
+                    nodecount += 1
+            if nodecount > 1:
+                refedge_fem_volumeelements.append(elem)
+        # for every refedge_fem_volumeelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
+        for elem in refedge_fem_volumeelements:
+            fe_refedge_nodes = []
+            for node in femelement_table[elem]:
+                if node in refedge_nodes:
+                    fe_refedge_nodes.append(node)
+                edge_table[elem] = fe_refedge_nodes  # { volumeID : ( edgenodeID, ... , edgenodeID  )} # only the refedge nodes
+        #  FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
+        edge_table = delete_duplicate_mesh_elements(edge_table)
+    elif is_face_femmesh(femmesh):
+        refedge_fem_faceelements = []
+        # if at least two nodes of a femfaceelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
+        for elem in femelement_table:
+            nodecount = 0
+            for node in femelement_table[elem]:
+                if node in refedge_nodes:
+                    nodecount += 1
+            if nodecount > 1:
+                refedge_fem_faceelements.append(elem)
+        # for every refedge_fem_faceelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
+        for elem in refedge_fem_faceelements:
+            fe_refedge_nodes = []
+            for node in femelement_table[elem]:
+                if node in refedge_nodes:
+                    fe_refedge_nodes.append(node)
+                edge_table[elem] = fe_refedge_nodes  # { faceID : ( edgenodeID, ... , edgenodeID  )} # only the refedge nodes
+        #  FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
+        edge_table = delete_duplicate_mesh_elements(edge_table)
+    elif is_edge_femmesh(femmesh):
+        refedge_fem_edgeelements = get_femelements_by_femnodes(femelement_table, refedge_nodes)
+        for elem in refedge_fem_edgeelements:
+            edge_table[elem] = femelement_table[elem]  # { edgeID : ( nodeID, ... , nodeID  )} # all nodes off this femedgeelement
+    return edge_table
+
+
+def get_ref_edgenodes_lengths(femnodes_mesh, edge_table):
+    # calulate the appropriate node_length for every node of every mesh edge (me)
+    # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
+
+    #  [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
+    if (not femnodes_mesh) or (not edge_table):
+        print('Error: empty femnodes_mesh or edge_table')
+        return []
+    node_length_table = []
+    mesh_edge_length = 0
+    # print(len(edge_table))
+    for me in edge_table:
+        femmesh_edgetype = len(edge_table[me])
+        if femmesh_edgetype == 2:  # 2 node femmesh edge
+            # end_node_length = mesh_edge_length / 2
+            #    ______
+            #  P1      P2
+            P1 = femnodes_mesh[edge_table[me][0]]
+            P2 = femnodes_mesh[edge_table[me][1]]
+            edge_vec = P2 - P1
+            mesh_edge_length = edge_vec.Length
+            # print(mesh_edge_length)
+            end_node_length = mesh_edge_length / 2.0
+            node_length_table.append((edge_table[me][0], end_node_length))
+            node_length_table.append((edge_table[me][1], end_node_length))
+
+        elif femmesh_edgetype == 3:  # 3 node femmesh edge
+            # end_node_length = mesh_edge_length / 6
+            # middle_node_length = mesh_face_area * 2 / 3
+            #   _______ _______
+            # P1       P3      P2
+            P1 = femnodes_mesh[edge_table[me][0]]
+            P2 = femnodes_mesh[edge_table[me][1]]
+            P3 = femnodes_mesh[edge_table[me][2]]
+            edge_vec1 = P3 - P1
+            edge_vec2 = P2 - P3
+            mesh_edge_length = edge_vec1.Length + edge_vec2.Length
+            # print(me, ' --> ', mesh_edge_length)
+            end_node_length = mesh_edge_length / 6.0
+            middle_node_length = mesh_edge_length * 2.0 / 3.0
+            node_length_table.append((edge_table[me][0], end_node_length))
+            node_length_table.append((edge_table[me][1], end_node_length))
+            node_length_table.append((edge_table[me][2], middle_node_length))
+    return node_length_table
+
+
+def get_ref_facenodes_table(femmesh, femelement_table, ref_face):
+    face_table = {}  # { meshfaceID : ( nodeID, ... , nodeID ) }
+    if is_solid_femmesh(femmesh):
+        if has_no_face_data(femmesh):
+            # there is no face data, the volumeID is used as key { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
+            ref_face_volume_elements = femmesh.getccxVolumesByFace(ref_face)  # list of tupels (mv, ccx_face_nr)
+            ref_face_nodes = femmesh.getNodesByFace(ref_face)
+            for ve in ref_face_volume_elements:
+                veID = ve[0]
+                ve_ref_face_nodes = []
+                for nodeID in femelement_table[veID]:
+                    if nodeID in ref_face_nodes:
+                        ve_ref_face_nodes.append(nodeID)
+                face_table[veID] = ve_ref_face_nodes  # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
+        else:  # the femmesh has face_data
+            volume_faces = femmesh.getVolumesByFace(ref_face)   # (mv, mf)
+            for mv, mf in volume_faces:
+                face_table[mf] = femmesh.getElementNodes(mf)
+    elif is_face_femmesh(femmesh):
+        ref_face_nodes = femmesh.getNodesByFace(ref_face)
+        ref_face_elements = get_femelements_by_femnodes(femelement_table, ref_face_nodes)
+        for mf in ref_face_elements:
+            face_table[mf] = femelement_table[mf]
+    return face_table
+
+
+def get_ref_facenodes_areas(femnodes_mesh, face_table):
+    # calulate the appropriate node_areas for every node of every mesh face (mf)
+    # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
+    # FIXME only gives exact results in case of a real triangle. If for S6 or C3D10 elements
+    # the midnodes are not on the line between the end nodes the area will not be a triangle
+    # see http://forum.freecadweb.org/viewtopic.php?f=18&t=10939&start=40#p91355  and ff
+
+    #  [ (nodeID,Area), ... , (nodeID,Area) ]  some nodes will have more than one entry
+    if (not femnodes_mesh) or (not face_table):
+        print('Error: empty femnodes_mesh or face_table')
+        return []
+    node_area_table = []
+    mesh_face_area = 0
+    for mf in face_table:
+        femmesh_facetype = len(face_table[mf])
+        if femmesh_facetype == 3:  # 3 node femmesh face triangle
+            # corner_node_area = mesh_face_area / 3.0
+            #      P3
+            #      /\
+            #     /  \
+            #    /____\
+            #  P1      P2
+            P1 = femnodes_mesh[face_table[mf][0]]
+            P2 = femnodes_mesh[face_table[mf][1]]
+            P3 = femnodes_mesh[face_table[mf][2]]
+
+            mesh_face_area = get_triangle_area(P1, P2, P3)
+            corner_node_area = mesh_face_area / 3.0
+
+            node_area_table.append((face_table[mf][0], corner_node_area))
+            node_area_table.append((face_table[mf][1], corner_node_area))
+            node_area_table.append((face_table[mf][2], corner_node_area))
+
+        elif femmesh_facetype == 4:  # 4 node femmesh face quad
+            FreeCAD.Console.PrintError('Face load on 4 node quad faces are not supported\n')
+
+        elif femmesh_facetype == 6:  # 6 node femmesh face triangle
+            # corner_node_area = 0
+            # middle_node_area = mesh_face_area / 3.0
+            #         P3
+            #         /\
+            #        /t3\
+            #       /    \
+            #     P6------P5
+            #     / \ t4 / \
+            #    /t1 \  /t2 \
+            #   /_____\/_____\
+            # P1      P4      P2
+            P1 = femnodes_mesh[face_table[mf][0]]
+            P2 = femnodes_mesh[face_table[mf][1]]
+            P3 = femnodes_mesh[face_table[mf][2]]
+            P4 = femnodes_mesh[face_table[mf][3]]
+            P5 = femnodes_mesh[face_table[mf][4]]
+            P6 = femnodes_mesh[face_table[mf][5]]
+
+            mesh_face_t1_area = get_triangle_area(P1, P4, P6)
+            mesh_face_t2_area = get_triangle_area(P2, P5, P4)
+            mesh_face_t3_area = get_triangle_area(P3, P6, P5)
+            mesh_face_t4_area = get_triangle_area(P4, P5, P6)
+            mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area
+            middle_node_area = mesh_face_area / 3.0
+
+            node_area_table.append((face_table[mf][0], 0))
+            node_area_table.append((face_table[mf][1], 0))
+            node_area_table.append((face_table[mf][2], 0))
+            node_area_table.append((face_table[mf][3], middle_node_area))
+            node_area_table.append((face_table[mf][4], middle_node_area))
+            node_area_table.append((face_table[mf][5], middle_node_area))
+
+        elif femmesh_facetype == 8:  # 8 node femmesh face quad
+            FreeCAD.Console.PrintError('Face load on 8 node quad faces are not supported\n')
+    return node_area_table
+
+
+def get_ref_shape_node_sum_geom_table(node_geom_table):
+    # shape could be Edge or Face, geom could be lenght or area
+    # summ of legth or area for each node of the ref_shape
+    node_sum_geom_table = {}
+    for n, A in node_geom_table:
+        # print(n, ' --> ', A)
+        if n in node_sum_geom_table:
+            node_sum_geom_table[n] = node_sum_geom_table[n] + A
+        else:
+            node_sum_geom_table[n] = A
+    return node_sum_geom_table
+
+
+def femelements_count_ok(femelement_table, count_femelements):
+    if count_femelements == len(femelement_table):
+        # print('Count FEM elements for the calculated node load distribution: ', count_femelements)
+        # print('Count FEM elements of the FreeCAD FEM mesh:  ', len(femelement_table))
+        return True
+    else:
+        print('ERROR: femelement_table != count_femelements')
+        print('Count FEM elements for the calculated node load distribution: ', count_femelements)
+        print('Count FEM Elements of the FreeCAD FEM Mesh:  ', len(femelement_table))
+        return False
+
+
+def delete_duplicate_mesh_elements(refelement_table):
+    new_refelement_table = {}  # duplicates deleted
+    for elem, nodes in refelement_table.items():
+        if sorted(nodes) not in sortlistoflistvalues(new_refelement_table.values()):
+            new_refelement_table[elem] = nodes
+    return new_refelement_table
+
+
+def get_triangle_area(P1, P2, P3):
+    vec1 = P2 - P1
+    vec2 = P3 - P1
+    vec3 = vec1.cross(vec2)
+    return 0.5 * vec3.Length
+
+
+def sortlistoflistvalues(listoflists):
+    new_list = []
+    for l in listoflists:
+        new_list.append(sorted(l))
+    return new_list
+
+
+def is_solid_femmesh(femmesh):
+    if femmesh.VolumeCount > 0:  # solid femmesh
+        return True
+
+
+def has_no_face_data(femmesh):
+    if femmesh.FaceCount == 0:   # femmesh has no face data, could be a edge femmesh or a solid femmesh without face data
+        return True
+
+
+def is_face_femmesh(femmesh):
+    if femmesh.VolumeCount == 0 and femmesh.FaceCount > 0:  # face femmesh
+        return True
+
+
+def is_edge_femmesh(femmesh):
+    if femmesh.VolumeCount == 0 and femmesh.FaceCount == 0 and femmesh.EdgeCount > 0:  # edge femmesh
+        return True
+
+
+def make_femmesh(mesh_data):
+    ''' makes an FreeCAD FEM Mesh object from FEM Mesh data
+    '''
+    import Fem
+    mesh = Fem.FemMesh()
+    m = mesh_data
+    if ('Nodes' in m) and (len(m['Nodes']) > 0):
+        if (('Hexa8Elem' in m) or ('Penta6Elem' in m) or ('Tetra4Elem' in m) or ('Tetra10Elem' in m) or
+           ('Penta6Elem' in m) or ('Hexa20Elem' in m) or ('Tria3Elem' in m) or ('Tria6Elem' in m) or
+           ('Quad4Elem' in m) or ('Quad8Elem' in m) or ('Seg2Elem' in m)):
+            nds = m['Nodes']
+
+            for i in nds:
+                n = nds[i]
+                mesh.addNode(n[0], n[1], n[2], i)
+            elms_hexa8 = m['Hexa8Elem']
+            for i in elms_hexa8:
+                e = elms_hexa8[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
+            elms_penta6 = m['Penta6Elem']
+            for i in elms_penta6:
+                e = elms_penta6[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5]], i)
+            elms_tetra4 = m['Tetra4Elem']
+            for i in elms_tetra4:
+                e = elms_tetra4[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3]], i)
+            elms_tetra10 = m['Tetra10Elem']
+            for i in elms_tetra10:
+                e = elms_tetra10[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9]], i)
+            elms_penta15 = m['Penta15Elem']
+            for i in elms_penta15:
+                e = elms_penta15[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
+                                e[10], e[11], e[12], e[13], e[14]], i)
+            elms_hexa20 = m['Hexa20Elem']
+            for i in elms_hexa20:
+                e = elms_hexa20[i]
+                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
+                                e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19]], i)
+            elms_tria3 = m['Tria3Elem']
+            for i in elms_tria3:
+                e = elms_tria3[i]
+                mesh.addFace([e[0], e[1], e[2]], i)
+            elms_tria6 = m['Tria6Elem']
+            for i in elms_tria6:
+                e = elms_tria6[i]
+                mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5]], i)
+            elms_quad4 = m['Quad4Elem']
+            for i in elms_quad4:
+                e = elms_quad4[i]
+                mesh.addFace([e[0], e[1], e[2], e[3]], i)
+            elms_quad8 = m['Quad8Elem']
+            for i in elms_quad8:
+                e = elms_quad8[i]
+                mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
+            elms_seg2 = m['Seg2Elem']
+            for i in elms_seg2:
+                e = elms_seg2[i]
+                mesh.addEdge(e[0], e[1])
+            print("imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15".format(
+                  len(nds), len(elms_hexa8), len(elms_penta6), len(elms_tetra4), len(elms_tetra10), len(elms_penta15)))
+            print("imported mesh: {} HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2".format(
+                  len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2)))
+    else:
+        FreeCAD.Console.PrintError("No Nodes found!\n")
+    return mesh

src/Mod/Fem/ccxFrdReader.py

@@ -300,7 +300,6 @@ def importFrd(filename, analysis=None):
     m = readResult(filename)
     mesh_object = None
     if(len(m['Nodes']) > 0):
-        import Fem
         if analysis is None:
             analysis_name = os.path.splitext(os.path.basename(filename))[0]
             import FemAnalysis
@@ -321,67 +320,11 @@ def importFrd(filename, analysis=None):
             z_span = abs(p_z_max - p_z_min)
             span = max(x_span, y_span, z_span)
 
-        if (not analysis) and ('Nodes' in m) and \
-            (('Hexa8Elem' in m) or ('Penta6Elem' in m) or ('Tetra4Elem' in m) or ('Tetra10Elem' in m) or
-             ('Penta6Elem' in m) or ('Hexa20Elem' in m) or ('Tria3Elem' in m) or ('Tria6Elem' in m) or
-             ('Quad4Elem' in m) or ('Quad8Elem' in m) or ('Seg2Elem' in m)):
-            mesh = Fem.FemMesh()
-            nds = m['Nodes']
+        if (not analysis):
+            import FemMeshTools
+            mesh = FemMeshTools.make_femmesh(m)
 
-            for i in nds:
-                n = nds[i]
-                mesh.addNode(n[0], n[1], n[2], i)
-            elms_hexa8 = m['Hexa8Elem']
-            for i in elms_hexa8:
-                e = elms_hexa8[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
-            elms_penta6 = m['Penta6Elem']
-            for i in elms_penta6:
-                e = elms_penta6[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5]], i)
-            elms_tetra4 = m['Tetra4Elem']
-            for i in elms_tetra4:
-                e = elms_tetra4[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3]], i)
-            elms_tetra10 = m['Tetra10Elem']
-            for i in elms_tetra10:
-                e = elms_tetra10[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9]], i)
-            elms_penta15 = m['Penta15Elem']
-            for i in elms_penta15:
-                e = elms_penta15[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
-                                e[10], e[11], e[12], e[13], e[14]], i)
-            elms_hexa20 = m['Hexa20Elem']
-            for i in elms_hexa20:
-                e = elms_hexa20[i]
-                mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
-                                e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19]], i)
-            elms_tria3 = m['Tria3Elem']
-            for i in elms_tria3:
-                e = elms_tria3[i]
-                mesh.addFace([e[0], e[1], e[2]], i)
-            elms_tria6 = m['Tria6Elem']
-            for i in elms_tria6:
-                e = elms_tria6[i]
-                mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5]], i)
-            elms_quad4 = m['Quad4Elem']
-            for i in elms_quad4:
-                e = elms_quad4[i]
-                mesh.addFace([e[0], e[1], e[2], e[3]], i)
-            elms_quad8 = m['Quad8Elem']
-            for i in elms_quad8:
-                e = elms_quad8[i]
-                mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
-            elms_seg2 = m['Seg2Elem']
-            for i in elms_seg2:
-                e = elms_seg2[i]
-                mesh.addEdge(e[0], e[1])
-            print ("imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15".format(
-                   len(nds), len(elms_hexa8), len(elms_penta6), len(elms_tetra4), len(elms_tetra10), len(elms_penta15)))
-            print ("imported mesh: {} HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2".format(
-                   len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2)))
-            if len(nds) > 0:
+            if len(m['Nodes']) > 0:
                 mesh_object = FreeCAD.ActiveDocument.addObject('Fem::FemMeshObject', 'ResultMesh')
                 mesh_object.FemMesh = mesh
                 analysis_object.Member = analysis_object.Member + [mesh_object]
@@ -432,8 +375,8 @@ def importFrd(filename, analysis=None):
                     results.StressValues = mstress
 
             if (results.NodeNumbers != 0 and results.NodeNumbers != stress.keys()):
-                print ("Inconsistent FEM results: element number for Stress doesn't equal element number for Displacement {} != {}"
-                       .format(results.NodeNumbers, len(results.StressValues)))
+                print("Inconsistent FEM results: element number for Stress doesn't equal element number for Displacement {} != {}"
+                      .format(results.NodeNumbers, len(results.StressValues)))
                 results.NodeNumbers = stress.keys()
 
             x_min, y_min, z_min = map(min, zip(*displacement))

src/Mod/Fem/ccxInpWriter.py

@@ -26,6 +26,7 @@ import FreeCAD
 import os
 import sys
 import time
+import FemMeshTools
 
 
 __title__ = "ccxInpWriter"
@@ -64,10 +65,12 @@ class inp_writer:
         self.fc_ver = FreeCAD.Version()
         self.ccx_eall = 'Eall'
         self.ccx_elsets = []
-        self.fem_mesh_nodes = {}
+        self.femmesh = self.mesh_object.FemMesh
+        self.femnodes_mesh = {}
+        self.femelement_table = {}
 
     def write_calculix_input_file(self):
-        self.mesh_object.FemMesh.writeABAQUS(self.file_name)
+        self.femmesh.writeABAQUS(self.file_name)
 
         # reopen file with "append" and add the analysis definition
         inpfile = open(self.file_name, 'a')
@@ -131,83 +134,78 @@ class inp_writer:
                 f.write('**No elements found for these objects\n')
 
     def write_node_sets_constraints_fixed(self, f):
+        # get nodes
+        for femobj in self.fixed_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
+        # write nodes to file
         f.write('\n***********************************************************\n')
         f.write('** Node set for fixed constraint\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for fobj in self.fixed_objects:
-            fix_obj = fobj['Object']
-            f.write('*NSET,NSET=' + fix_obj.Name + '\n')
-            for o, elem in fix_obj.References:
-                fo = o.Shape.getElement(elem)
-                n = []
-                if fo.ShapeType == 'Face':
-                    n = self.mesh_object.FemMesh.getNodesByFace(fo)
-                elif fo.ShapeType == 'Edge':
-                    n = self.mesh_object.FemMesh.getNodesByEdge(fo)
-                elif fo.ShapeType == 'Vertex':
-                    n = self.mesh_object.FemMesh.getNodesByVertex(fo)
-                for i in n:
-                    f.write(str(i) + ',\n')
+        for femobj in self.fixed_objects:
+            f.write('*NSET,NSET=' + femobj['Object'].Name + '\n')
+            for n in femobj['Nodes']:
+                f.write(str(n) + ',\n')
 
     def write_node_sets_constraints_displacement(self, f):
+        # get nodes
+        for femobj in self.displacement_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
+        # write nodes to file
         f.write('\n***********************************************************\n')
         f.write('** Node sets for prescribed displacement constraint\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for fobj in self.displacement_objects:
-            disp_obj = fobj['Object']
-            f.write('*NSET,NSET=' + disp_obj.Name + '\n')
-            for o, elem in disp_obj.References:
-                fo = o.Shape.getElement(elem)
-                n = []
-                if fo.ShapeType == 'Face':
-                    n = self.mesh_object.FemMesh.getNodesByFace(fo)
-                elif fo.ShapeType == 'Edge':
-                    n = self.mesh_object.FemMesh.getNodesByEdge(fo)
-                elif fo.ShapeType == 'Vertex':
-                    n = self.mesh_object.FemMesh.getNodesByVertex(fo)
-                for i in n:
-                    f.write(str(i) + ',\n')
+        for femobj in self.displacement_objects:
+            f.write('*NSET,NSET=' + femobj['Object'].Name + '\n')
+            for n in femobj['Nodes']:
+                f.write(str(n) + ',\n')
 
     def write_node_sets_constraints_force(self, f):
-        f.write('\n***********************************************************\n')
-        f.write('** Node sets for loads\n')
-        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for fobj in self.force_objects:
-            frc_obj = fobj['Object']
+        # check shape type of reference shape
+        for femobj in self.force_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            frc_obj = femobj['Object']
             # in GUI defined frc_obj all ref_shape have the same shape type
             # TODO in FemTools: check if all RefShapes really have the same type an write type to dictionary
-            fobj['RefShapeType'] = ''
+            femobj['RefShapeType'] = ''
             if frc_obj.References:
                 first_ref_obj = frc_obj.References[0]
                 first_ref_shape = first_ref_obj[0].Shape.getElement(first_ref_obj[1])
-                fobj['RefShapeType'] = first_ref_shape.ShapeType
+                femobj['RefShapeType'] = first_ref_shape.ShapeType
             else:
                 # frc_obj.References could be empty ! # TODO in FemTools: check
                 FreeCAD.Console.PrintError('At least one Force Object has empty References!\n')
-            if fobj['RefShapeType'] == 'Vertex':
-                pass  # point load on vertices --> the FemElementTable is not needed for node load calculation
-            elif fobj['RefShapeType'] == 'Face' and is_solid_mesh(self.mesh_object.FemMesh) and not has_no_face_data(self.mesh_object.FemMesh):
-                pass  # solid_mesh with face data --> the FemElementTable is not needed for node load calculation
+            if femobj['RefShapeType'] == 'Vertex':
+                print("load on vertices --> we do not need the femelement_table and femnodes_mesh for node load calculation")
+                pass
+            elif femobj['RefShapeType'] == 'Face' and FemMeshTools.is_solid_mesh(self.femmesh) and not FemMeshTools.has_no_face_data(self.femmesh):
+                print("solid_mesh with face data --> we do not need the femelement_table but we need the femnodes_mesh for node load calculation")
+                if not self.femnodes_mesh:
+                    self.femnodes_mesh = self.femmesh.Nodes
             else:
-                if not hasattr(self, 'fem_element_table'):
-                    self.fem_element_table = getFemElementTable(self.mesh_object.FemMesh)
-        for fobj in self.force_objects:
-            if fobj['RefShapeType'] == 'Vertex':
-                frc_obj = fobj['Object']
-                f.write('*NSET,NSET=' + frc_obj.Name + '\n')
-                NbrForceNodes = 0
-                for o, elem in frc_obj.References:
-                    fo = o.Shape.getElement(elem)
-                    n = []
-                    n = self.mesh_object.FemMesh.getNodesByVertex(fo)
-                    for i in n:
-                        f.write(str(i) + ',\n')
-                        NbrForceNodes = NbrForceNodes + 1   # NodeSum of mesh-nodes of ALL reference shapes from force_object
+                print("mesh without needed data --> we need the femelement_table and femnodes_mesh for node load calculation")
+                if not self.femnodes_mesh:
+                    self.femnodes_mesh = self.femmesh.Nodes
+                if not self.femelement_table:
+                    self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
+        # point loads on vertices, get nodes
+        for femobj in self.force_objects:
+            if femobj['RefShapeType'] == 'Vertex':
+                femobj['Nodes'] = FemMeshTools.get_femnodes_by_references(self.femmesh, femobj['Object'].References)
                 # calculate node load
-                if NbrForceNodes != 0:
-                    fobj['NodeLoad'] = (frc_obj.Force) / NbrForceNodes
-                    f.write('** concentrated load [N] distributed on all mesh nodes of the given vertieces\n')
-                    f.write('** ' + str(frc_obj.Force) + ' N / ' + str(NbrForceNodes) + ' Nodes = ' + str(fobj['NodeLoad']) + ' N on each node\n')
+                if len(femobj['Nodes']) != 0:
+                    femobj['NodeLoad'] = (femobj['Object'].Force) / len(femobj['Nodes'])
+        # point loads on vertices, write nodes to file
+        f.write('\n***********************************************************\n')
+        f.write('** Node sets for loads\n')
+        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
+        for femobj in self.force_objects:
+            if femobj['RefShapeType'] == 'Vertex':
+                frc_obj = femobj['Object']
+                f.write('*NSET,NSET=' + frc_obj.Name + '\n')
+                for n in femobj['Nodes']:
+                    f.write(str(n) + ',\n')
+                f.write('** concentrated load [N] distributed on all mesh nodes of the given vertieces\n')
+                f.write('** ' + str(frc_obj.Force) + ' N / ' + str(len(femobj['Nodes'])) + ' Nodes = ' + str(femobj['NodeLoad']) + ' N on each node\n')
+                f.write('\n')
             else:
                 f.write('** no point load on vertices --> no set for node loads\n')
 
@@ -216,8 +214,8 @@ class inp_writer:
         f.write('** Materials\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
         f.write('** Young\'s modulus unit is MPa = N/mm2\n')
-        for m in self.material_objects:
-            mat_obj = m['Object']
+        for femobj in self.material_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            mat_obj = femobj['Object']
             # get material properties
             YM = FreeCAD.Units.Quantity(mat_obj.Material['YoungsModulus'])
             YM_in_MPa = YM.getValueAs('MPa')
@@ -275,8 +273,8 @@ class inp_writer:
         f.write('\n***********************************************************\n')
         f.write('** Constaints\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for fixed_object in self.fixed_objects:
-            fix_obj_name = fixed_object['Object'].Name
+        for femobj in self.fixed_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            fix_obj_name = femobj['Object'].Name
             f.write('*BOUNDARY\n')
             f.write(fix_obj_name + ',1\n')
             f.write(fix_obj_name + ',2\n')
@@ -291,52 +289,59 @@ class inp_writer:
         f.write('\n***********************************************************\n')
         f.write('** Displacement constraint applied\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for disp_obj in self.displacement_objects:
-            disp_obj_name = disp_obj['Object'].Name
+        for femobj in self.displacement_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            disp_obj = femobj['Object']
+            disp_obj_name = disp_obj.Name
             f.write('*BOUNDARY\n')
-            if disp_obj['Object'].xFix:
+            if disp_obj.xFix:
                 f.write(disp_obj_name + ',1\n')
-            elif not disp_obj['Object'].xFree:
-                f.write(disp_obj_name + ',1,1,' + str(disp_obj['Object'].xDisplacement) + '\n')
-            if disp_obj['Object'].yFix:
+            elif not disp_obj.xFree:
+                f.write(disp_obj_name + ',1,1,' + str(disp_obj.xDisplacement) + '\n')
+            if disp_obj.yFix:
                 f.write(disp_obj_name + ',2\n')
-            elif not disp_obj['Object'].yFree:
-                f.write(disp_obj_name + ',2,2,' + str(disp_obj['Object'].yDisplacement) + '\n')
-            if disp_obj['Object'].zFix:
+            elif not disp_obj.yFree:
+                f.write(disp_obj_name + ',2,2,' + str(disp_obj.yDisplacement) + '\n')
+            if disp_obj.zFix:
                 f.write(disp_obj_name + ',3\n')
-            elif not disp_obj['Object'].zFree:
-                f.write(disp_obj_name + ',3,3,' + str(disp_obj['Object'].zDisplacement) + '\n')
+            elif not disp_obj.zFree:
+                f.write(disp_obj_name + ',3,3,' + str(disp_obj.zDisplacement) + '\n')
 
             if self.beamsection_objects or self.shellthickness_objects:
-                if disp_obj['Object'].rotxFix:
+                if disp_obj.rotxFix:
                     f.write(disp_obj_name + ',4\n')
-                elif not disp_obj['Object'].rotxFree:
-                    f.write(disp_obj_name + ',4,4,' + str(disp_obj['Object'].xRotation) + '\n')
-                if disp_obj['Object'].rotyFix:
+                elif not disp_obj.rotxFree:
+                    f.write(disp_obj_name + ',4,4,' + str(disp_obj.xRotation) + '\n')
+                if disp_obj.rotyFix:
                     f.write(disp_obj_name + ',5\n')
-                elif not disp_obj['Object'].rotyFree:
-                    f.write(disp_obj_name + ',5,5,' + str(disp_obj['Object'].yRotation) + '\n')
-                if disp_obj['Object'].rotzFix:
+                elif not disp_obj.rotyFree:
+                    f.write(disp_obj_name + ',5,5,' + str(disp_obj.yRotation) + '\n')
+                if disp_obj.rotzFix:
                     f.write(disp_obj_name + ',6\n')
-                elif not disp_obj['Object'].rotzFree:
-                    f.write(disp_obj_name + ',6,6,' + str(disp_obj['Object'].zRotation) + '\n')
+                elif not disp_obj.rotzFree:
+                    f.write(disp_obj_name + ',6,6,' + str(disp_obj.zRotation) + '\n')
         f.write('\n')
 
     def write_constraints_force(self, f):
+        # get node loads
+        for femobj in self.force_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            frc_obj = femobj['Object']
+            if frc_obj.Force == 0:
+                print('  Warning --> Force = 0')
+            if femobj['RefShapeType'] == 'Edge':  # line load on edges
+                femobj['NodeLoadTable'] = FemMeshTools.get_force_obj_edge_nodeload_table(self.femmesh, self.femelement_table, self.femnodes_mesh, frc_obj)
+            elif femobj['RefShapeType'] == 'Face':  # area load on faces
+                femobj['NodeLoadTable'] = FemMeshTools.get_force_obj_face_nodeload_table(self.femmesh, self.femelement_table, self.femnodes_mesh, frc_obj)
+        # write node loads to file
         f.write('\n***********************************************************\n')
         f.write('** Node loads\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
         f.write('*CLOAD\n')
-        for fobj in self.force_objects:
-            frc_obj = fobj['Object']
-            f.write('** ' + frc_obj.Name + '\n')
-            direction_vec = frc_obj.DirectionVector
-            if frc_obj.Force == 0:
-                print('  Warning --> Force = 0')
-
-            if fobj['RefShapeType'] == 'Vertex':  # point load on vertieces
-                node_load = fobj['NodeLoad']
-                frc_obj_name = frc_obj.Name
+        for femobj in self.force_objects:  # femobj --> dict, FreeCAD document object is femobj['Object']
+            frc_obj_name = femobj['Object'].Name
+            direction_vec = femobj['Object'].DirectionVector
+            f.write('** ' + frc_obj_name + '\n')
+            if femobj['RefShapeType'] == 'Vertex':  # point load on vertieces
+                node_load = femobj['NodeLoad']
                 f.write('** force: ' + str(node_load) + ' N,  direction: ' + str(direction_vec) + '\n')
                 v1 = "{:.13E}".format(direction_vec.x * node_load)
                 v2 = "{:.13E}".format(direction_vec.y * node_load)
@@ -344,47 +349,11 @@ class inp_writer:
                 f.write(frc_obj_name + ',1,' + v1 + '\n')
                 f.write(frc_obj_name + ',2,' + v2 + '\n')
                 f.write(frc_obj_name + ',3,' + v3 + '\n\n')
-
-            elif fobj['RefShapeType'] == 'Edge':  # line load on edges
-                sum_ref_edge_length = 0
-                sum_ref_edge_node_length = 0  # for debugging
-                sum_node_load = 0  # for debugging
-                for o, elem in frc_obj.References:
-                    elem_o = o.Shape.getElement(elem)
-                    sum_ref_edge_length += elem_o.Length
-                if sum_ref_edge_length != 0:
-                    force_per_sum_ref_edge_length = frc_obj.Force / sum_ref_edge_length
-                for o, elem in frc_obj.References:
-                    elem_o = o.Shape.getElement(elem)
-                    ref_edge = elem_o
-
-                    # edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
-                    edge_table = self.get_refedge_node_table(ref_edge)
-
-                    # node_length_table = [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
-                    node_length_table = self.get_refedge_node_lengths(edge_table)
-
-                    # node_sum_length_table = { nodeID : Length, ... , nodeID : Length }  LengthSum for each node, one entry for each node
-                    node_sum_length_table = self.get_ref_shape_node_sum_geom_table(node_length_table)
-
-                    # node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad }  NodeLoad for each node, one entry for each node
-                    node_load_table = {}
-                    sum_node_lengths = 0  # for debugging
-                    for node in node_sum_length_table:
-                        sum_node_lengths += node_sum_length_table[node]  # for debugging
-                        node_load_table[node] = node_sum_length_table[node] * force_per_sum_ref_edge_length
-                    ratio_refedge_lengths = sum_node_lengths / elem_o.Length
-                    if ratio_refedge_lengths < 0.99 or ratio_refedge_lengths > 1.01:
-                        FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
-                        print('  sum_node_lengths:', sum_node_lengths)
-                        print('  refedge_length:  ', elem_o.Length)
-                        bad_refedge = elem_o
-                    sum_ref_edge_node_length += sum_node_lengths
-
-                    f.write('** node loads on element ' + fobj['RefShapeType'] + ': ' + o.Name + ':' + elem + '\n')
-                    for n in sorted(node_load_table):
-                        node_load = node_load_table[n]
-                        sum_node_load += node_load  # for debugging
+            elif femobj['RefShapeType'] == 'Edge' or femobj['RefShapeType'] == 'Face':
+                for ref_shape in femobj['NodeLoadTable']:
+                    f.write('** ' + ref_shape[0] + '\n')
+                    for n in sorted(ref_shape[1]):
+                        node_load = ref_shape[1][n]
                         if (direction_vec.x != 0.0):
                             v1 = "{:.13E}".format(direction_vec.x * node_load)
                             f.write(str(n) + ',1,' + v1 + '\n')
@@ -396,125 +365,19 @@ class inp_writer:
                             f.write(str(n) + ',3,' + v3 + '\n')
                     f.write('\n')
                 f.write('\n')
-                ratio = sum_node_load / frc_obj.Force
-                if ratio < 0.99 or ratio > 1.01:
-                    print('Deviation  sum_node_load to frc_obj.Force is more than 1% :  ', ratio)
-                    print('  sum_ref_edge_node_length: ', sum_ref_edge_node_length)
-                    print('  sum_ref_edge_length:      ', sum_ref_edge_length)
-                    print('  sum_node_load:          ', sum_node_load)
-                    print('  frc_obj.Force:          ', frc_obj.Force)
-                    print('  the reason could be simply a circle length --> see method get_ref_edge_node_lengths')
-                    print('  the reason could also be an problem in retrieving the ref_edge_node_length')
-
-                    # try debugging of the last bad refedge
-                    print('DEBUGGING')
-                    print(bad_refedge)
-
-                    print('bad_refedge_nodes')
-                    bad_refedge_nodes = self.mesh_object.FemMesh.getNodesByEdge(bad_refedge)
-                    print(len(bad_refedge_nodes))
-                    print(bad_refedge_nodes)
-                    # import FreeCADGui
-                    # FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_refedge_nodes
-
-                    print('bad_edge_table')
-                    # bad_edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
-                    bad_edge_table = self.get_refedge_node_table(bad_refedge)
-                    print(len(bad_edge_table))
-                    bad_edge_table_nodes = []
-                    for elem in bad_edge_table:
-                        print(elem, ' --> ', bad_edge_table[elem])
-                        for node in bad_edge_table[elem]:
-                            if node not in bad_edge_table_nodes:
-                                bad_edge_table_nodes.append(node)
-                    print('sorted(bad_edge_table_nodes)')
-                    print(sorted(bad_edge_table_nodes))   # should be == bad_refedge_nodes
-                    # import FreeCADGui
-                    # FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_edge_table_nodes
-                    # bad_node_length_table = [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
-
-                    print('good_edge_table')
-                    good_edge_table = delete_duplicate_mesh_elements(bad_edge_table)
-                    for elem in good_edge_table:
-                        print(elem, ' --> ', bad_edge_table[elem])
-
-                    print('bad_node_length_table')
-                    bad_node_length_table = self.get_refedge_node_lengths(bad_edge_table)
-                    for n, l in bad_node_length_table:
-                        print(n, ' --> ', l)
-
-            elif fobj['RefShapeType'] == 'Face':  # area load on faces
-                sum_ref_face_area = 0
-                sum_ref_face_node_area = 0  # for debugging
-                sum_node_load = 0  # for debugging
-                for o, elem in frc_obj.References:
-                    elem_o = o.Shape.getElement(elem)
-                    sum_ref_face_area += elem_o.Area
-                if sum_ref_face_area != 0:
-                    force_per_sum_ref_face_area = frc_obj.Force / sum_ref_face_area
-                for o, elem in frc_obj.References:
-                    elem_o = o.Shape.getElement(elem)
-                    ref_face = elem_o
-
-                    # face_table = { meshfaceID : ( nodeID, ... , nodeID ) }
-                    face_table = self.get_ref_face_node_table(ref_face)
-
-                    # node_area_table = [ (nodeID, Area), ... , (nodeID, Area) ]  some nodes will have more than one entry
-                    node_area_table = self.get_ref_face_node_areas(face_table)
-
-                    # node_sum_area_table = { nodeID : Area, ... , nodeID : Area }  AreaSum for each node, one entry for each node
-                    node_sum_area_table = self.get_ref_shape_node_sum_geom_table(node_area_table)
-
-                    # node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad }  NodeLoad for each node, one entry for each node
-                    node_load_table = {}
-                    sum_node_areas = 0  # for debugging
-                    for node in node_sum_area_table:
-                        sum_node_areas += node_sum_area_table[node]  # for debugging
-                        node_load_table[node] = node_sum_area_table[node] * force_per_sum_ref_face_area
-                    ratio_refface_areas = sum_node_areas / elem_o.Area
-                    if ratio_refface_areas < 0.99 or ratio_refface_areas > 1.01:
-                        FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
-                        print('  sum_node_lengths:', sum_node_areas)
-                        print('  refedge_length:  ', elem_o.Area)
-                    sum_ref_face_node_area += sum_node_areas
-
-                    f.write('** node loads on element ' + fobj['RefShapeType'] + ': ' + o.Name + ':' + elem + '\n')
-                    for n in sorted(node_load_table):
-                        node_load = node_load_table[n]
-                        sum_node_load += node_load  # for debugging
-                        if (direction_vec.x != 0.0):
-                            v1 = "{:.13E}".format(direction_vec.x * node_load)
-                            f.write(str(n) + ',1,' + v1 + '\n')
-                        if (direction_vec.y != 0.0):
-                            v2 = "{:.13E}".format(direction_vec.y * node_load)
-                            f.write(str(n) + ',2,' + v2 + '\n')
-                        if (direction_vec.z != 0.0):
-                            v3 = "{:.13E}".format(direction_vec.z * node_load)
-                            f.write(str(n) + ',3,' + v3 + '\n')
-                    f.write('\n')
-                f.write('\n')
-                ratio = sum_node_load / frc_obj.Force
-                if ratio < 0.99 or ratio > 1.01:
-                    print('Deviation  sum_node_load to frc_obj.Force is more than 1% :  ', ratio)
-                    print('  sum_ref_face_node_area: ', sum_ref_face_node_area)
-                    print('  sum_ref_face_area:      ', sum_ref_face_area)
-                    print('  sum_node_load:          ', sum_node_load)
-                    print('  frc_obj.Force:          ', frc_obj.Force)
-                    print('  the reason could be simply a circle area --> see method get_ref_face_node_areas')
-                    print('  the reason could also be an problem in retrieving the ref_face_node_area')
 
     def write_constraints_pressure(self, f):
         f.write('\n***********************************************************\n')
         f.write('** Element + CalculiX face + load in [MPa]\n')
         f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
-        for fobj in self.pressure_objects:
-            prs_obj = fobj['Object']
+        for femobj in self.pressure_objects:
+            prs_obj = femobj['Object']
             f.write('*DLOAD\n')
             for o, e in prs_obj.References:
                 rev = -1 if prs_obj.Reversed else 1
                 elem = o.Shape.getElement(e)
                 if elem.ShapeType == 'Face':
-                    v = self.mesh_object.FemMesh.getccxVolumesByFace(elem)
+                    v = self.femmesh.getccxVolumesByFace(elem)
                     f.write("** Load on face {}\n".format(e))
                     for i in v:
                         f.write("{},P{},{}\n".format(i[0], i[1], rev * prs_obj.Pressure))
@@ -607,8 +470,10 @@ class inp_writer:
         self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_single_mat_multiple_beam(self):
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
         mat_obj = self.material_objects[0]['Object']
-        self.get_femelement_sets(self.beamsection_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.beamsection_objects)
         for beamsec_data in self.beamsection_objects:
             beamsec_obj = beamsec_data['Object']
             ccx_elset = {}
@@ -620,8 +485,10 @@ class inp_writer:
             self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_single_mat_multiple_shell(self):
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
         mat_obj = self.material_objects[0]['Object']
-        self.get_femelement_sets(self.shellthickness_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.shellthickness_objects)
         for shellth_data in self.shellthickness_objects:
             shellth_obj = shellth_data['Object']
             ccx_elset = {}
@@ -633,8 +500,10 @@ class inp_writer:
             self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_multiple_mat_single_beam(self):
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
         beamsec_obj = self.beamsection_objects[0]['Object']
-        self.get_femelement_sets(self.material_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.material_objects)
         for mat_data in self.material_objects:
             mat_obj = mat_data['Object']
             ccx_elset = {}
@@ -646,8 +515,10 @@ class inp_writer:
             self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_multiple_mat_single_shell(self):
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
         shellth_obj = self.shellthickness_objects[0]['Object']
-        self.get_femelement_sets(self.material_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.material_objects)
         for mat_data in self.material_objects:
             mat_obj = mat_data['Object']
             ccx_elset = {}
@@ -659,7 +530,9 @@ class inp_writer:
             self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_multiple_mat_solid(self):
-        self.get_femelement_sets(self.material_objects)
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.material_objects)
         for mat_data in self.material_objects:
             mat_obj = mat_data['Object']
             ccx_elset = {}
@@ -670,8 +543,10 @@ class inp_writer:
             self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_multiple_mat_multiple_beam(self):
-        self.get_femelement_sets(self.beamsection_objects)
-        self.get_femelement_sets(self.material_objects)
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.beamsection_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.material_objects)
         for beamsec_data in self.beamsection_objects:
             beamsec_obj = beamsec_data['Object']
             for mat_data in self.material_objects:
@@ -689,8 +564,10 @@ class inp_writer:
                 self.ccx_elsets.append(ccx_elset)
 
     def get_ccx_elsets_multiple_mat_multiple_shell(self):
-        self.get_femelement_sets(self.shellthickness_objects)
-        self.get_femelement_sets(self.material_objects)
+        if not self.femelement_table:
+            self.femelement_table = FemMeshTools.get_femelement_table(self.femmesh)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.shellthickness_objects)
+        FemMeshTools.get_femelement_sets(self.femmesh, self.femelement_table, self.material_objects)
         for shellth_data in self.shellthickness_objects:
             shellth_obj = shellth_data['Object']
             for mat_data in self.material_objects:
@@ -707,344 +584,8 @@ class inp_writer:
                 ccx_elset['ccx_mat_name'] = mat_obj.Material['Name']
                 self.ccx_elsets.append(ccx_elset)
 
-    def get_femelement_sets(self, fem_objects):
-        # get femelements for reference shapes of each obj.References
-        if not hasattr(self, 'fem_element_table'):
-            self.fem_element_table = getFemElementTable(self.mesh_object.FemMesh)
-        count_femelements = 0
-        referenced_femelements = []
-        has_remaining_femelements = None
-        for fem_object_i, fem_object in enumerate(fem_objects):
-            obj = fem_object['Object']
-            fem_object['ShortName'] = get_ccx_elset_short_name(obj, fem_object_i)  # unique short ccx_identifier
-            if obj.References:
-                ref_shape_femelements = []
-                for ref in obj.References:
-                    femnodes = []
-                    femelements = []
-                    if ref[1]:
-                        r = ref[0].Shape.getElement(ref[1])
-                    else:
-                        r = ref[0].Shape
-                    # print('  ReferenceShape : ', r.ShapeType, ', ', ref[0].Name, ', ', ref[0].Label, ' --> ', ref[1])
-                    if r.ShapeType == 'Edge':
-                        femnodes = self.mesh_object.FemMesh.getNodesByEdge(r)
-                    elif r.ShapeType == 'Face':
-                        femnodes = self.mesh_object.FemMesh.getNodesByFace(r)
-                    elif r.ShapeType == 'Solid':
-                        femnodes = self.mesh_object.FemMesh.getNodesBySolid(r)
-                    else:
-                        print('  No Edge, Face or Solid as reference shapes!')
-                    femelements = getFemElementsByNodes(self.fem_element_table, femnodes)
-                    ref_shape_femelements += femelements
-                    referenced_femelements += femelements
-                    count_femelements += len(femelements)
-                fem_object['FEMElements'] = ref_shape_femelements
-            else:
-                has_remaining_femelements = obj.Name
-        # get remaining femelements for the fem_objects
-        if has_remaining_femelements:
-            remaining_femelements = []
-            for elemid in self.fem_element_table:
-                if elemid not in referenced_femelements:
-                    remaining_femelements.append(elemid)
-            count_femelements += len(remaining_femelements)
-            for fem_object in fem_objects:
-                obj = fem_object['Object']
-                if obj.Name == has_remaining_femelements:
-                    fem_object['FEMElements'] = sorted(remaining_femelements)
-        # check if all worked out well
-        if not femelements_count_ok(self.fem_element_table, count_femelements):
-            FreeCAD.Console.PrintError('Error in get_femelement_sets -- > femelements_count_ok failed!\n')
-
-    def get_refedge_node_table(self, refedge):
-        edge_table = {}  # { meshedgeID : ( nodeID, ... , nodeID ) }
-        refedge_nodes = self.mesh_object.FemMesh.getNodesByEdge(refedge)
-        if is_solid_mesh(self.mesh_object.FemMesh):
-            refedge_fem_volumeelements = []
-            # if at least two nodes of a femvolumeelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
-            for elem in self.fem_element_table:
-                nodecount = 0
-                for node in self.fem_element_table[elem]:
-                    if node in refedge_nodes:
-                        nodecount += 1
-                if nodecount > 1:
-                    refedge_fem_volumeelements.append(elem)
-            # for every refedge_fem_volumeelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
-            for elem in refedge_fem_volumeelements:
-                fe_refedge_nodes = []
-                for node in self.fem_element_table[elem]:
-                    if node in refedge_nodes:
-                        fe_refedge_nodes.append(node)
-                    edge_table[elem] = fe_refedge_nodes  # { volumeID : ( edgenodeID, ... , edgenodeID  )} # only the refedge nodes
-            #  FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
-            edge_table = delete_duplicate_mesh_elements(edge_table)
-        elif is_shell_mesh(self.mesh_object.FemMesh):
-            refedge_fem_faceelements = []
-            # if at least two nodes of a femfaceelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
-            for elem in self.fem_element_table:
-                nodecount = 0
-                for node in self.fem_element_table[elem]:
-                    if node in refedge_nodes:
-                        nodecount += 1
-                if nodecount > 1:
-                    refedge_fem_faceelements.append(elem)
-            # for every refedge_fem_faceelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
-            for elem in refedge_fem_faceelements:
-                fe_refedge_nodes = []
-                for node in self.fem_element_table[elem]:
-                    if node in refedge_nodes:
-                        fe_refedge_nodes.append(node)
-                    edge_table[elem] = fe_refedge_nodes  # { faceID : ( edgenodeID, ... , edgenodeID  )} # only the refedge nodes
-            #  FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
-            edge_table = delete_duplicate_mesh_elements(edge_table)
-        elif is_beam_mesh(self.mesh_object.FemMesh):
-            refedge_fem_edgeelements = getFemElementsByNodes(self.fem_element_table, refedge_nodes)
-            for elem in refedge_fem_edgeelements:
-                edge_table[elem] = self.fem_element_table[elem]  # { edgeID : ( nodeID, ... , nodeID  )} # all nodes off this femedgeelement
-        return edge_table
-
-    def get_ref_face_node_table(self, ref_face):
-        face_table = {}  # { meshfaceID : ( nodeID, ... , nodeID ) }
-        if is_solid_mesh(self.mesh_object.FemMesh):
-            if has_no_face_data(self.mesh_object.FemMesh):
-                # there is no face data, the volumeID is used as key { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
-                ref_face_volume_elements = self.mesh_object.FemMesh.getccxVolumesByFace(ref_face)  # list of tupels (mv, ccx_face_nr)
-                ref_face_nodes = self.mesh_object.FemMesh.getNodesByFace(ref_face)
-                for ve in ref_face_volume_elements:
-                    veID = ve[0]
-                    ve_ref_face_nodes = []
-                    for nodeID in self.fem_element_table[veID]:
-                        if nodeID in ref_face_nodes:
-                            ve_ref_face_nodes.append(nodeID)
-                    face_table[veID] = ve_ref_face_nodes  # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
-            else:  # the femmesh has face_data
-                volume_faces = self.mesh_object.FemMesh.getVolumesByFace(ref_face)   # (mv, mf)
-                for mv, mf in volume_faces:
-                    face_table[mf] = self.mesh_object.FemMesh.getElementNodes(mf)
-        elif is_shell_mesh(self.mesh_object.FemMesh):
-            ref_face_nodes = self.mesh_object.FemMesh.getNodesByFace(ref_face)
-            ref_face_elements = getFemElementsByNodes(self.fem_element_table, ref_face_nodes)
-            for mf in ref_face_elements:
-                face_table[mf] = self.fem_element_table[mf]
-        return face_table
-
-    def get_refedge_node_lengths(self, edge_table):
-        # calulate the appropriate node_length for every node of every mesh edge (me)
-        # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
-
-        #  [ (nodeID, length), ... , (nodeID, length) ]  some nodes will have more than one entry
-        node_length_table = []
-        mesh_edge_length = 0
-        if not self.fem_mesh_nodes:
-            self.fem_mesh_nodes = self.mesh_object.FemMesh.Nodes
-        # print(len(edge_table))
-        for me in edge_table:
-            if len(edge_table[me]) == 2:  # 2 node mesh edge
-                # end_node_length = mesh_edge_length / 2
-                #    ______
-                #  P1      P2
-                P1 = self.fem_mesh_nodes[edge_table[me][0]]
-                P2 = self.fem_mesh_nodes[edge_table[me][1]]
-                edge_vec = P2 - P1
-                mesh_edge_length = edge_vec.Length
-                # print(mesh_edge_length)
-                end_node_length = mesh_edge_length / 2.0
-                node_length_table.append((edge_table[me][0], end_node_length))
-                node_length_table.append((edge_table[me][1], end_node_length))
-
-            elif len(edge_table[me]) == 3:  # 3 node mesh edge
-                # end_node_length = mesh_edge_length / 6
-                # middle_node_length = mesh_face_area * 2 / 3
-                #   _______ _______
-                # P1       P3      P2
-                P1 = self.fem_mesh_nodes[edge_table[me][0]]
-                P2 = self.fem_mesh_nodes[edge_table[me][1]]
-                P3 = self.fem_mesh_nodes[edge_table[me][2]]
-                edge_vec1 = P3 - P1
-                edge_vec2 = P2 - P3
-                mesh_edge_length = edge_vec1.Length + edge_vec2.Length
-                # print(me, ' --> ', mesh_edge_length)
-                end_node_length = mesh_edge_length / 6.0
-                middle_node_length = mesh_edge_length * 2.0 / 3.0
-                node_length_table.append((edge_table[me][0], end_node_length))
-                node_length_table.append((edge_table[me][1], end_node_length))
-                node_length_table.append((edge_table[me][2], middle_node_length))
-        return node_length_table
-
-    def get_ref_face_node_areas(self, face_table):
-        # calulate the appropriate node_areas for every node of every mesh face (mf)
-        # G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
-        # FIXME only gives exact results in case of a real triangle. If for S6 or C3D10 elements
-        # the midnodes are not on the line between the end nodes the area will not be a triangle
-        # see http://forum.freecadweb.org/viewtopic.php?f=18&t=10939&start=40#p91355  and ff
-
-        #  [ (nodeID,Area), ... , (nodeID,Area) ]  some nodes will have more than one entry
-        node_area_table = []
-        mesh_face_area = 0
-        if not self.fem_mesh_nodes:
-            self.fem_mesh_nodes = self.mesh_object.FemMesh.Nodes
-        for mf in face_table:
-            if len(face_table[mf]) == 3:  # 3 node mesh face triangle
-                # corner_node_area = mesh_face_area / 3.0
-                #      P3
-                #      /\
-                #     /  \
-                #    /____\
-                #  P1      P2
-                P1 = self.fem_mesh_nodes[face_table[mf][0]]
-                P2 = self.fem_mesh_nodes[face_table[mf][1]]
-                P3 = self.fem_mesh_nodes[face_table[mf][2]]
-
-                mesh_face_area = getTriangleArea(P1, P2, P3)
-                corner_node_area = mesh_face_area / 3.0
-
-                node_area_table.append((face_table[mf][0], corner_node_area))
-                node_area_table.append((face_table[mf][1], corner_node_area))
-                node_area_table.append((face_table[mf][2], corner_node_area))
-
-            elif len(face_table[mf]) == 4:  # 4 node mesh face quad
-                FreeCAD.Console.PrintError('Face load on 4 node quad faces are not supported\n')
-
-            elif len(face_table[mf]) == 6:  # 6 node mesh face triangle
-                # corner_node_area = 0
-                # middle_node_area = mesh_face_area / 3.0
-                #         P3
-                #         /\
-                #        /t3\
-                #       /    \
-                #     P6------P5
-                #     / \ t4 / \
-                #    /t1 \  /t2 \
-                #   /_____\/_____\
-                # P1      P4      P2
-                P1 = self.fem_mesh_nodes[face_table[mf][0]]
-                P2 = self.fem_mesh_nodes[face_table[mf][1]]
-                P3 = self.fem_mesh_nodes[face_table[mf][2]]
-                P4 = self.fem_mesh_nodes[face_table[mf][3]]
-                P5 = self.fem_mesh_nodes[face_table[mf][4]]
-                P6 = self.fem_mesh_nodes[face_table[mf][5]]
-
-                mesh_face_t1_area = getTriangleArea(P1, P4, P6)
-                mesh_face_t2_area = getTriangleArea(P2, P5, P4)
-                mesh_face_t3_area = getTriangleArea(P3, P6, P5)
-                mesh_face_t4_area = getTriangleArea(P4, P5, P6)
-                mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area
-                middle_node_area = mesh_face_area / 3.0
-
-                node_area_table.append((face_table[mf][0], 0))
-                node_area_table.append((face_table[mf][1], 0))
-                node_area_table.append((face_table[mf][2], 0))
-                node_area_table.append((face_table[mf][3], middle_node_area))
-                node_area_table.append((face_table[mf][4], middle_node_area))
-                node_area_table.append((face_table[mf][5], middle_node_area))
-
-            elif len(face_table[mf]) == 8:  # 8 node mesh face quad
-                FreeCAD.Console.PrintError('Face load on 8 node quad faces are not supported\n')
-        return node_area_table
-
-    def get_ref_shape_node_sum_geom_table(self, node_geom_table):
-        # shape could be Edge or Face, geom could be lenght or area
-        # summ of legth or area for each node of the ref_shape
-        node_sum_geom_table = {}
-        for n, A in node_geom_table:
-            # print(n, ' --> ', A)
-            if n in node_sum_geom_table:
-                node_sum_geom_table[n] = node_sum_geom_table[n] + A
-            else:
-                node_sum_geom_table[n] = A
-        return node_sum_geom_table
-
 
 # Helpers
-def getTriangleArea(P1, P2, P3):
-    vec1 = P2 - P1
-    vec2 = P3 - P1
-    vec3 = vec1.cross(vec2)
-    return 0.5 * vec3.Length
-
-
-def getFemElementTable(fem_mesh):
-    """ getFemElementTable(fem_mesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }"""
-    fem_element_table = {}
-    if is_solid_mesh(fem_mesh):
-        for i in fem_mesh.Volumes:
-            fem_element_table[i] = fem_mesh.getElementNodes(i)
-    elif is_shell_mesh(fem_mesh):
-        for i in fem_mesh.Faces:
-            fem_element_table[i] = fem_mesh.getElementNodes(i)
-    elif is_beam_mesh(fem_mesh):
-        for i in fem_mesh.Edges:
-            fem_element_table[i] = fem_mesh.getElementNodes(i)
-    else:
-        FreeCAD.Console.PrintError('Neither solid nor shell nor beam mesh!\n')
-    return fem_element_table
-
-
-def getFemElementsByNodes(fem_element_table, node_list):
-    '''for every fem_element of fem_element_table
-    if all nodes of the fem_element are in node_list,
-    the fem_element is added to the list which is returned
-    e: elementlist
-    nodes: nodelist '''
-    e = []  # elementlist
-    for elementID in sorted(fem_element_table):
-        nodecount = 0
-        for nodeID in fem_element_table[elementID]:
-            if nodeID in node_list:
-                nodecount = nodecount + 1
-        if nodecount == len(fem_element_table[elementID]):   # all nodes of the element are in the node_list!
-            e.append(elementID)
-    return e
-
-
-def is_solid_mesh(fem_mesh):
-    if fem_mesh.VolumeCount > 0:  # solid mesh
-        return True
-
-
-def has_no_face_data(fem_mesh):
-    if fem_mesh.FaceCount == 0:   # mesh has no face data, could be a beam mesh or a solid mesh without face data
-        return True
-
-
-def is_shell_mesh(fem_mesh):
-    if fem_mesh.VolumeCount == 0 and fem_mesh.FaceCount > 0:  # shell mesh
-        return True
-
-
-def is_beam_mesh(fem_mesh):
-    if fem_mesh.VolumeCount == 0 and fem_mesh.FaceCount == 0 and fem_mesh.EdgeCount > 0:  # beam mesh
-        return True
-
-
-def femelements_count_ok(fem_element_table, count_femelements):
-    if count_femelements == len(fem_element_table):
-        # print('Count Elements written to CalculiX file: ', count_femelements)
-        # print('Count Elements of the FreeCAD FEM Mesh:  ', len(fem_element_table))
-        return True
-    else:
-        print('ERROR: self.fem_element_table != count_femelements')
-        print('Count Elements written to CalculiX file: ', count_femelements)
-        print('Count Elements of the FreeCAD FEM Mesh:  ', len(fem_element_table))
-        return False
-
-
-def delete_duplicate_mesh_elements(refelement_table):
-    new_refelement_table = {}  # duplicates deleted
-    for elem, nodes in refelement_table.items():
-        if sorted(nodes) not in sortlistoflistvalues(new_refelement_table.values()):
-            new_refelement_table[elem] = nodes
-    return new_refelement_table
-
-
-def sortlistoflistvalues(listoflists):
-    new_list = []
-    for l in listoflists:
-        new_list.append(sorted(l))
-    return new_list
-
-
 def get_ccx_elset_beam_name(mat_name, beamsec_name, mat_short_name=None, beamsec_short_name=None):
     if not mat_short_name:
         mat_short_name = 'Mat0'
@@ -1076,14 +617,3 @@ def get_ccx_elset_solid_name(mat_name, solid_name=None, mat_short_name=None):
         return mat_short_name + solid_name
     else:
         return mat_name + solid_name
-
-
-def get_ccx_elset_short_name(obj, i):
-    if hasattr(obj, "Proxy") and obj.Proxy.Type == 'MechanicalMaterial':
-        return 'Mat' + str(i)
-    elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemBeamSection':
-        return 'Beam' + str(i)
-    elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemShellThickness':
-        return 'Shell' + str(i)
-    else:
-        print('Error: ', obj.Name, ' --> ', obj.Proxy.Type)

src/Mod/Fem/test_files/ccx/cube_static.inp

@@ -468,11 +468,6 @@ Eall
 ** Node sets for prescribed displacement constraint
 ** written by write_node_sets_constraints_displacement function
 
-***********************************************************
-** Node sets for loads
-** written by write_node_sets_constraints_force function
-** no point load on vertices --> no set for node loads
-
 ***********************************************************
 ** Materials
 ** written by write_materials function
@@ -515,7 +510,7 @@ FemConstraintFixed,3
 ** written by write_constraints_force function
 *CLOAD
 ** FemConstraintForce
-** node loads on element Face: Box:Face6
+** node loads on shape: Box:Face6
 2,3,-0.0000000000000E+00
 4,3,-0.0000000000000E+00
 6,3,-0.0000000000000E+00