FEM: Add frequency analysis to FEM wb

Signed-off-by: Przemo Firszt <przemo@firszt.eu>

src/Mod/Fem/FemTools.py

@@ -33,6 +33,8 @@ class FemTools(QtCore.QRunnable, QtCore.QObject):
         QtCore.QRunnable.__init__(self)
         QtCore.QObject.__init__(self)
 
+        self.known_analysis_types = ["static", "frequency"]
+
         if analysis:
             self.analysis = analysis
         else:
@@ -40,6 +42,8 @@ class FemTools(QtCore.QRunnable, QtCore.QObject):
             self.analysis = FemGui.getActiveAnalysis()
         if self.analysis:
             self.update_objects()
+            self.set_analysis_type()
+            self.set_eigenmode_parameters()
             self.base_name = ""
             self.results_present = False
             self.setup_working_dir()
@@ -119,12 +123,15 @@ class FemTools(QtCore.QRunnable, QtCore.QObject):
         message = ""
         if not self.analysis:
             message += "No active Analysis\n"
+        if self.analysis_type not in self.known_analysis_types:
+            message += "Unknown analysis type: {}\n".format(self.analysis_type)
         if not self.mesh:
             message += "No mesh object in the Analysis\n"
         if not self.material:
             message += "No material object in the Analysis\n"
         if not self.fixed_constraints:
             message += "No fixed-constraint nodes defined in the Analysis\n"
+        if self.analysis_type == "static":
             if not (self.force_constraints or self.pressure_constraints):
                 message += "No force-constraint or pressure-constraint defined in the Analysis\n"
         return message
@@ -136,7 +143,8 @@ class FemTools(QtCore.QRunnable, QtCore.QObject):
         try:
             inp_writer = iw.inp_writer(self.analysis, self.mesh, self.material,
                                        self.fixed_constraints, self.force_constraints,
-                                       self.pressure_constraints, self.working_dir)
+                                       self.pressure_constraints, self.analysis_type,
+                                       self.eigenmode_parameters, self.working_dir)
             self.base_name = inp_writer.write_calculix_input_file()
         except:
             print "Unexpected error when writing CalculiX input file:", sys.exc_info()[0]
@@ -165,6 +173,26 @@ class FemTools(QtCore.QRunnable, QtCore.QObject):
             return p.returncode
         return -1
 
+    ## sets eigenmode parameters for CalculiX frequency analysis
+    #  @param self The python object self
+    #  @param number number of eigenmodes that wll be calculated, default 10
+    #  @param limit_low lower value of requested eigenfrequency range, default 0.0
+    #  @param limit_high higher value of requested eigenfrequency range, default 1000000.0
+    def set_eigenmode_parameters(self, number=10, limit_low=0.0, limit_high=1000000.0):
+        self.eigenmode_parameters = (number, limit_low, limit_high)
+
+    def set_base_name(self, base_name=None):
+        if base_name is None:
+            self.base_name = ""
+        else:
+            self.base_name = base_name
+
+    def set_analysis_type(self, analysis_type=None):
+        if analysis_type is None:
+            self.analysis_type = "static"
+        else:
+            self.analysis_type = analysis_type
+
     def setup_working_dir(self, working_dir=None):
         if working_dir is None:
             self.fem_prefs = FreeCAD.ParamGet("User parameter:BaseApp/Preferences/Mod/Fem")

src/Mod/Fem/Gui/Resources/Fem.qrc

@@ -17,6 +17,7 @@
         <file>icons/fem-new-analysis.svg</file>
         <file>icons/fem-purge-results.svg</file>
         <file>icons/fem-quick-analysis.svg</file>
+        <file>icons/fem-frequency-analysis.svg</file>
         <file>icons/fem-result.svg</file>
         <file>icons/preferences-fem.svg</file>
         <file>translations/Fem_af.qm</file>

src/Mod/Fem/Gui/Workbench.cpp

@@ -70,6 +70,7 @@ Gui::ToolBarItem* Workbench::setupToolBars() const
          << "Separator"
          << "Fem_MechanicalJobControl"
          << "Fem_Quick_Analysis"
+         << "Fem_Frequency_Analysis"
          << "Fem_PurgeResults"
          << "Fem_ShowResult";
     return root;
@@ -97,6 +98,7 @@ Gui::MenuItem* Workbench::setupMenuBar() const
          << "Separator"
          << "Fem_MechanicalJobControl"
          << "Fem_Quick_Analysis"
+         << "Fem_Frequency_Analysis"
          << "Fem_PurgeResults"
          << "Fem_ShowResult";
 

src/Mod/Fem/MechanicalAnalysis.py

@@ -179,6 +179,36 @@ class _CommandQuickAnalysis:
         return FreeCADGui.ActiveDocument is not None and FemGui.getActiveAnalysis() is not None
 
 
+class _CommandFrequencyAnalysis:
+    def GetResources(self):
+        return {'Pixmap': 'fem-frequency-analysis',
+                'MenuText': QtCore.QT_TRANSLATE_NOOP("Fem_Frequency_Analysis", "Run frequency analysis with CalculiX ccx"),
+                'Accel': "R, F",
+                'ToolTip': QtCore.QT_TRANSLATE_NOOP("Fem_Frequency_Analysis", "Write .inp file and run frequency analysis with CalculiX ccx")}
+
+    def Activated(self):
+        def load_results(ret_code):
+            if ret_code == 0:
+                self.fea.load_results()
+            else:
+                print "CalculiX failed ccx finished with error {}".format(ret_code)
+
+        self.fea = FemTools()
+        self.fea.purge_results()
+        self.fea.reset_mesh_color()
+        self.fea.reset_mesh_deformation()
+        self.fea.set_analysis_type('frequency')
+        message = self.fea.check_prerequisites()
+        if message:
+            QtGui.QMessageBox.critical(None, "Missing prerequisite", message)
+            return
+        self.fea.finished.connect(load_results)
+        QtCore.QThreadPool.globalInstance().start(self.fea)
+
+    def IsActive(self):
+        return FreeCADGui.ActiveDocument is not None and FemGui.getActiveAnalysis() is not None
+
+
 class _CommandMechanicalShowResult:
     "the Fem JobControl command definition"
     def GetResources(self):
@@ -595,19 +625,28 @@ class _ResultControlTaskPanel:
         self.form.le_max.setProperty("unit", unit)
         self.form.le_max.setText("{:.6} {}".format(maxm, unit))
 
+    def update_displacement(self, factor=None):
+        if factor is None:
+            if FreeCAD.FEM_dialog["show_disp"]:
+                factor = self.form.hsb_displacement_factor.value()
+            else:
+                factor = 0
+        self.MeshObject.ViewObject.applyDisplacement(factor)
+
     def show_displacement(self, checked):
         QApplication.setOverrideCursor(Qt.WaitCursor)
         FreeCAD.FEM_dialog["show_disp"] = checked
-        factor = 0.0
+        if "result_object" in FreeCAD.FEM_dialog:
-        if checked:
+            if FreeCAD.FEM_dialog["result_object"] != self.result_object:
-            factor = self.form.hsb_displacement_factor.value()
+                self.update_displacement(reset=True)
+        FreeCAD.FEM_dialog["result_object"] = self.result_object
         self.MeshObject.ViewObject.setNodeDisplacementByVectors(self.result_object.ElementNumbers, self.result_object.DisplacementVectors)
-        self.MeshObject.ViewObject.applyDisplacement(factor)
+        self.update_displacement()
         QtGui.qApp.restoreOverrideCursor()
 
     def hsb_disp_factor_changed(self, value):
-        self.MeshObject.ViewObject.applyDisplacement(value)
         self.form.sb_displacement_factor.setValue(value)
+        self.update_displacement()
 
     def sb_disp_factor_max_changed(self, value):
         FreeCAD.FEM_dialog["disp_factor_max"] = value
@@ -659,5 +698,6 @@ if FreeCAD.GuiUp:
     FreeCADGui.addCommand('Fem_CreateFromShape', _CommandFemFromShape())
     FreeCADGui.addCommand('Fem_MechanicalJobControl', _CommandMechanicalJobControl())
     FreeCADGui.addCommand('Fem_Quick_Analysis', _CommandQuickAnalysis())
+    FreeCADGui.addCommand('Fem_Frequency_Analysis', _CommandFrequencyAnalysis())
     FreeCADGui.addCommand('Fem_PurgeResults', _CommandPurgeFemResults())
     FreeCADGui.addCommand('Fem_ShowResult', _CommandMechanicalShowResult())

src/Mod/Fem/ccxFrdReader.py

@@ -34,21 +34,21 @@ if open.__module__ == '__builtin__':
     pyopen = open  # because we'll redefine open below
 
 
-displacement = []
-
-
 # read a calculix result file and extract the nodes, displacement vectores and stress values.
 def readResult(frd_input):
     frd_file = pyopen(frd_input, "r")
     nodes = {}
-    disp = {}
-    stress = {}
     elements = {}
+    results = []
+    mode_results = {}
+    mode_disp = {}
+    mode_stress = {}
 
-    disp_found = False
+    mode_disp_found = False
     nodes_found = False
-    stress_found = False
+    mode_stress_found = False
     elements_found = False
+    eigenmode = 0
     elem = -1
     elemType = 0
 
@@ -83,21 +83,23 @@ def readResult(frd_input):
             node_id_8 = int(line[83:93])
             node_id_10 = int(line[93:103])
             elements[elem] = (node_id_1, node_id_2, node_id_3, node_id_4, node_id_5, node_id_6, node_id_7, node_id_8, node_id_9, node_id_10)
+        #Check if we found new eigenmode
+        if line[5:10] == "PMODE":
+            eigenmode = int(line[30:36])
         #Check if we found displacement section
         if line[5:9] == "DISP":
-            disp_found = True
+            mode_disp_found = True
         #we found a displacement line in the frd file
-        if disp_found and (line[1:3] == "-1"):
+        if mode_disp_found and (line[1:3] == "-1"):
             elem = int(line[4:13])
-            disp_x = float(line[13:25])
+            mode_disp_x = float(line[13:25])
-            disp_y = float(line[25:37])
+            mode_disp_y = float(line[25:37])
-            disp_z = float(line[37:49])
+            mode_disp_z = float(line[37:49])
-            disp[elem] = FreeCAD.Vector(disp_x, disp_y, disp_z)
+            mode_disp[elem] = FreeCAD.Vector(mode_disp_x, mode_disp_y, mode_disp_z)
-            displacement.append((disp_x, disp_y, disp_z))
         if line[5:11] == "STRESS":
-            stress_found = True
+            mode_stress_found = True
         #we found a displacement line in the frd file
-        if stress_found and (line[1:3] == "-1"):
+        if mode_stress_found and (line[1:3] == "-1"):
             elem = int(line[4:13])
             stress_1 = float(line[13:25])
             stress_2 = float(line[25:37])
@@ -105,26 +107,50 @@ def readResult(frd_input):
             stress_4 = float(line[49:61])
             stress_5 = float(line[61:73])
             stress_6 = float(line[73:85])
-            stress[elem] = (stress_1, stress_2, stress_3, stress_4, stress_5, stress_6)
+            mode_stress[elem] = (stress_1, stress_2, stress_3, stress_4, stress_5, stress_6)
         #Check for the end of a section
         if line[1:3] == "-3":
-            #the section with the displacements and the nodes ended
+            if mode_disp_found:
-            disp_found = False
+                mode_disp_found = False
+
+            if mode_stress_found:
+                mode_stress_found = False
+
+            if mode_disp and mode_stress:
+                mode_results = {}
+                mode_results['number'] = eigenmode
+                mode_results['disp'] = mode_disp
+                mode_results['stress'] = mode_stress
+                results.append(mode_results)
+                mode_disp = {}
+                mode_stress = {}
+                eigenmode = 0
             nodes_found = False
-            stress_found = False
             elements_found = False
 
     frd_file.close()
-    FreeCAD.Console.PrintLog('Read CalculiX result: {} Nodes, {} Displacements and {} Stress values\n'.format(len(nodes), len(disp), len(stress)))
+    return {'Nodes': nodes, 'Tet10Elem': elements, 'Results': results}
 
-    return {'Nodes': nodes, 'Tet10Elem': elements, 'Displacement': disp, 'Stress': stress}
+
+def calculate_von_mises(i):
+    # Von mises stress (http://en.wikipedia.org/wiki/Von_Mises_yield_criterion)
+    s11 = i[0]
+    s22 = i[1]
+    s33 = i[2]
+    s12 = i[3]
+    s23 = i[4]
+    s31 = i[5]
+    s11s22 = pow(s11 - s22, 2)
+    s22s33 = pow(s22 - s33, 2)
+    s33s11 = pow(s33 - s11, 2)
+    s12s23s31 = 6 * (pow(s12, 2) + pow(s23, 2) * pow(s31, 2))
+    vm_stress = sqrt(0.5 * (s11s22 + s22s33 + s33s11 + s12s23s31))
+    return vm_stress
 
 
 def importFrd(filename, Analysis=None):
-    mstress = []
-    global displacement
-    displacement = []
     m = readResult(filename)
+    result_set_number = len(m['Results'])
     MeshObject = None
     if(len(m) > 0):
         import Fem
@@ -134,11 +160,23 @@ def importFrd(filename, Analysis=None):
             AnalysisObject.Label = AnalysisName
         else:
             AnalysisObject = Analysis
-        results = FreeCAD.ActiveDocument.addObject('Fem::FemResultObject', 'Results')
+
+        if 'Nodes' in m:
+            positions = []
+            for k, v in m['Nodes'].iteritems():
+                positions.append(v)
+            p_x_max, p_y_max, p_z_max = map(max, zip(*positions))
+            p_x_min, p_y_min, p_z_min = map(min, zip(*positions))
+
+            x_span = abs(p_x_max - p_x_min)
+            y_span = abs(p_y_max - p_y_min)
+            z_span = abs(p_z_max - p_z_min)
+            span = max(x_span, y_span, z_span)
 
         if ('Tet10Elem' in m) and ('Nodes' in m) and (not Analysis):
             mesh = Fem.FemMesh()
             nds = m['Nodes']
+
             for i in nds:
                 n = nds[i]
                 mesh.addNode(n[0], n[1], n[2], i)
@@ -151,52 +189,68 @@ def importFrd(filename, Analysis=None):
                 MeshObject.FemMesh = mesh
                 AnalysisObject.Member = AnalysisObject.Member + [MeshObject]
 
-        if 'Displacement' in m:
+        for result_set in m['Results']:
-            disp = m['Displacement']
+            eigenmode_number = result_set['number']
+            if result_set_number > 1:
+                results_name = 'Mode_' + str(eigenmode_number) + '_results'
+            else:
+                results_name = 'Results'
+            results = FreeCAD.ActiveDocument.addObject('Fem::FemResultObject', results_name)
+
+            disp = result_set['disp']
+            l = len(disp)
+            displacement = []
+            for k, v in disp.iteritems():
+                displacement.append(v)
+
+            x_max, y_max, z_max = map(max, zip(*displacement))
+            if result_set_number > 1:
+                max_disp = max(x_max, y_max, z_max)
+                # Allow for max displacement to be 0.1% of the span
+                # FIXME - add to Preferences
+                max_allowed_disp = 0.001 * span
+                scale = max_allowed_disp / max_disp
+            else:
+                scale = 1.0
+
             if len(disp) > 0:
-                results.DisplacementVectors = disp.values()
+                results.DisplacementVectors = map((lambda x: x * scale), disp.values())
                 results.ElementNumbers = disp.keys()
                 if(MeshObject):
                     results.Mesh = MeshObject
-        if 'Stress' in m:
+
-            stress = m['Stress']
+            stress = result_set['stress']
             if len(stress) > 0:
+                mstress = []
                 for i in stress.values():
-                    # Von mises stress (http://en.wikipedia.org/wiki/Von_Mises_yield_criterion)
+                    mstress.append(calculate_von_mises(i))
-                    s11 = i[0]
+                if result_set_number > 1:
-                    s22 = i[1]
+                    results.StressValues = map((lambda x: x * scale), mstress)
-                    s33 = i[2]
+                else:
-                    s12 = i[3]
-                    s23 = i[4]
-                    s31 = i[5]
-                    s11s22 = pow(s11 - s22, 2)
-                    s22s33 = pow(s22 - s33, 2)
-                    s33s11 = pow(s33 - s11, 2)
-                    s12s23s31 = 6 * (pow(s12, 2) + pow(s23, 2) * pow(s31, 2))
-                    mstress.append(sqrt(0.5 * (s11s22 + s22s33 + s33s11 + s12s23s31)))
-
                     results.StressValues = mstress
-                if (results.ElementNumbers != 0 and results.ElementNumbers != stress.keys()):
-                    print "Inconsistent FEM results: element number for Stress doesn't equal element number for Displacement"
-                    results.ElementNumbers = stress.keys()
-                if(MeshObject):
-                    results.Mesh = MeshObject
 
-        l = len(displacement)
+            if (results.ElementNumbers != 0 and results.ElementNumbers != stress.keys()):
-        x_max, y_max, z_max = map(max, zip(*displacement))
+                print ("Inconsistent FEM results: element number for Stress doesn't equal element number for Displacement {} != {}"
+                       .format(results.ElementNumbers, len(results.StressValues)))
+                results.ElementNumbers = stress.keys()
+
             x_min, y_min, z_min = map(min, zip(*displacement))
             sum_list = map(sum, zip(*displacement))
             x_avg, y_avg, z_avg = [i / l for i in sum_list]
-        s_max = max(mstress)
+
-        s_min = min(mstress)
+            s_max = max(results.StressValues)
-        s_avg = sum(mstress) / l
+            s_min = min(results.StressValues)
+            s_avg = sum(results.StressValues) / l
+
             disp_abs = []
             for d in displacement:
                 disp_abs.append(sqrt(pow(d[0], 2) + pow(d[1], 2) + pow(d[2], 2)))
             results.DisplacementLengths = disp_abs
+
             a_max = max(disp_abs)
             a_min = min(disp_abs)
             a_avg = sum(disp_abs) / l
+
             results.Stats = [x_min, x_avg, x_max,
                              y_min, y_avg, y_max,
                              z_min, z_avg, z_max,

src/Mod/Fem/ccxInpWriter.py

@@ -5,7 +5,8 @@ import time
 
 
 class inp_writer:
-    def __init__(self, analysis_obj, mesh_obj, mat_obj, fixed_obj, force_obj, pressure_obj, dir_name=None):
+    def __init__(self, analysis_obj, mesh_obj, mat_obj, fixed_obj, force_obj,
+                 pressure_obj, analysis_type=None, eigenmode_parameters=None, dir_name=None):
         self.dir_name = dir_name
         self.analysis = analysis_obj
         self.mesh_object = mesh_obj
@@ -13,6 +14,11 @@ class inp_writer:
         self.fixed_objects = fixed_obj
         self.force_objects = force_obj
         self.pressure_objects = pressure_obj
+        if eigenmode_parameters:
+            self.no_of_eigenfrequencies = eigenmode_parameters[0]
+            self.eigenfrequeny_range_low = eigenmode_parameters[1]
+            self.eigenfrequeny_range_high = eigenmode_parameters[2]
+        self.analysis_type = analysis_type
         if not dir_name:
             self.dir_name = FreeCAD.ActiveDocument.TransientDir.replace('\\', '/') + '/FemAnl_' + analysis_obj.Uid[-4:]
         if not os.path.isdir(self.dir_name):
@@ -33,8 +39,11 @@ class inp_writer:
         self.write_materials(inpfile)
         self.write_step_begin(inpfile)
         self.write_constraints_fixed(inpfile)
+        if self.analysis_type is None or self.analysis_type == "static":
             self.write_constraints_force(inpfile)
             self.write_face_load(inpfile)
+        elif self.analysis_type == "frequency":
+            self.write_frequency(inpfile)
         self.write_outputs_types(inpfile)
         self.write_step_end(inpfile)
         self.write_footer(inpfile)
@@ -326,6 +335,13 @@ class inp_writer:
                     for i in v:
                         f.write("{},P{},{}\n".format(i[0], i[1], rev * prs_obj.Pressure))
 
+    def write_frequency(self, f):
+        f.write('\n***********************************************************\n')
+        f.write('** Frequency analysis\n')
+        f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
+        f.write('*FREQUENCY\n')
+        f.write('{},{},{}\n'.format(self.no_of_eigenfrequencies, self.eigenfrequeny_range_low, self.eigenfrequeny_range_high))
+
     def write_outputs_types(self, f):
         f.write('\n***********************************************************\n')
         f.write('** Outputs --> frd file\n')