diff --git a/src/Mod/Fem/FemTools.py b/src/Mod/Fem/FemTools.py index 4f768a7dfd..10cf5ad221 100644 --- a/src/Mod/Fem/FemTools.py +++ b/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,14 +123,17 @@ 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 not (self.force_constraints or self.pressure_constraints): - message += "No force-constraint or pressure-constraint 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 def write_inp_file(self): @@ -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") diff --git a/src/Mod/Fem/Gui/Resources/Fem.qrc b/src/Mod/Fem/Gui/Resources/Fem.qrc index 3d4c2ef3e2..9e25b250cc 100755 --- a/src/Mod/Fem/Gui/Resources/Fem.qrc +++ b/src/Mod/Fem/Gui/Resources/Fem.qrc @@ -17,6 +17,7 @@ icons/fem-new-analysis.svg icons/fem-purge-results.svg icons/fem-quick-analysis.svg + icons/fem-frequency-analysis.svg icons/fem-result.svg icons/preferences-fem.svg translations/Fem_af.qm diff --git a/src/Mod/Fem/Gui/Resources/icons/fem-frequency-analysis.svg b/src/Mod/Fem/Gui/Resources/icons/fem-frequency-analysis.svg new file mode 100644 index 0000000000..a92186aab1 --- /dev/null +++ b/src/Mod/Fem/Gui/Resources/icons/fem-frequency-analysis.svg @@ -0,0 +1,488 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + diff --git a/src/Mod/Fem/Gui/Workbench.cpp b/src/Mod/Fem/Gui/Workbench.cpp index af64fad6a6..0ca33fdbc3 100755 --- a/src/Mod/Fem/Gui/Workbench.cpp +++ b/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"; diff --git a/src/Mod/Fem/MechanicalAnalysis.py b/src/Mod/Fem/MechanicalAnalysis.py index bc72fe33a9..c68b7bd8dd 100644 --- a/src/Mod/Fem/MechanicalAnalysis.py +++ b/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 checked: - factor = self.form.hsb_displacement_factor.value() + if "result_object" in FreeCAD.FEM_dialog: + if FreeCAD.FEM_dialog["result_object"] != self.result_object: + 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()) diff --git a/src/Mod/Fem/ccxFrdReader.py b/src/Mod/Fem/ccxFrdReader.py index c79c9a80c3..4b1e37a55a 100644 --- a/src/Mod/Fem/ccxFrdReader.py +++ b/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]) - disp_y = float(line[25:37]) - disp_z = float(line[37:49]) - disp[elem] = FreeCAD.Vector(disp_x, disp_y, disp_z) - displacement.append((disp_x, disp_y, disp_z)) + mode_disp_x = float(line[13:25]) + mode_disp_y = float(line[25:37]) + mode_disp_z = float(line[37:49]) + mode_disp[elem] = FreeCAD.Vector(mode_disp_x, mode_disp_y, mode_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 - disp_found = False + if mode_disp_found: + 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,58 +189,74 @@ def importFrd(filename, Analysis=None): MeshObject.FemMesh = mesh AnalysisObject.Member = AnalysisObject.Member + [MeshObject] - if 'Displacement' in m: - disp = m['Displacement'] + for result_set in m['Results']: + 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) - 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)) - mstress.append(sqrt(0.5 * (s11s22 + s22s33 + s33s11 + s12s23s31))) + mstress.append(calculate_von_mises(i)) + if result_set_number > 1: + results.StressValues = map((lambda x: x * scale), mstress) + else: + results.StressValues = mstress - 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 + if (results.ElementNumbers != 0 and results.ElementNumbers != stress.keys()): + 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() - l = len(displacement) - x_max, y_max, z_max = map(max, zip(*displacement)) - 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_avg = sum(mstress) / 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, - a_min, a_avg, a_max, - s_min, s_avg, s_max] - AnalysisObject.Member = AnalysisObject.Member + [results] + 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(results.StressValues) + 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, + a_min, a_avg, a_max, + s_min, s_avg, s_max] + AnalysisObject.Member = AnalysisObject.Member + [results] if(FreeCAD.GuiUp): import FemGui diff --git a/src/Mod/Fem/ccxInpWriter.py b/src/Mod/Fem/ccxInpWriter.py index e80f41b794..fd97563722 100644 --- a/src/Mod/Fem/ccxInpWriter.py +++ b/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) - self.write_constraints_force(inpfile) - self.write_face_load(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')