diff --git a/src/Mod/Fem/importCcxFrdResults.py b/src/Mod/Fem/importCcxFrdResults.py index 17d3dc719..36b0710ff 100644 --- a/src/Mod/Fem/importCcxFrdResults.py +++ b/src/Mod/Fem/importCcxFrdResults.py @@ -32,8 +32,6 @@ __url__ = "http://www.freecadweb.org" import FreeCAD import os -from math import pow, sqrt -import numpy as np ########## generic FreeCAD import and export methods ########## @@ -63,6 +61,7 @@ def insert(filename, docname): ########## module specific methods ########## def importFrd(filename, analysis=None, result_name_prefix=None): + import importToolsFem import ObjectsFem if result_name_prefix is None: result_name_prefix = '' @@ -89,7 +88,6 @@ def importFrd(filename, analysis=None, result_name_prefix=None): span = max(x_span, y_span, z_span) if (not analysis): - import importToolsFem mesh = importToolsFem.make_femmesh(m) if len(m['Nodes']) > 0: @@ -108,171 +106,13 @@ def importFrd(filename, analysis=None, result_name_prefix=None): results_name = result_name_prefix + 'time_' + str(step_time) + '_results' else: results_name = result_name_prefix + 'results' + results = ObjectsFem.makeResultMechanical(results_name) - for m in analysis_object.Member: + for m in analysis_object.Member: # TODO analysis could have multiple mesh objects in the future if m.isDerivedFrom("Fem::FemMeshObject"): results.Mesh = m break - - try: - disp = result_set['disp'] - stressv = result_set['stressv'] - strainv = result_set['strainv'] - no_of_values = len(disp) - displacement = [] - for k, v in disp.items(): - displacement.append(v) - - x_max, y_max, z_max = map(max, zip(*displacement)) - if eigenmode_number > 0: - 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 = list(map((lambda x: x * scale), disp.values())) - results.StressVectors = list(map((lambda x: x * scale), stressv.values())) - results.StrainVectors = list(map((lambda x: x * scale), strainv.values())) - results.NodeNumbers = disp.keys() - if(mesh_object): - results.Mesh = mesh_object - - stress = result_set['stress'] - if len(stress) > 0: - mstress = [] - prinstress1 = [] - prinstress2 = [] - prinstress3 = [] - shearstress = [] - for i in stress.values(): - mstress.append(calculate_von_mises(i)) - prin1, prin2, prin3, shear = calculate_principal_stress(i) - prinstress1.append(prin1) - prinstress2.append(prin2) - prinstress3.append(prin3) - shearstress.append(shear) - if eigenmode_number > 0: - results.StressValues = list(map((lambda x: x * scale), mstress)) - results.PrincipalMax = list(map((lambda x: x * scale), prinstress1)) - results.PrincipalMed = list(map((lambda x: x * scale), prinstress2)) - results.PrincipalMin = list(map((lambda x: x * scale), prinstress3)) - results.MaxShear = list(map((lambda x: x * scale), shearstress)) - results.Eigenmode = eigenmode_number - else: - results.StressValues = mstress - results.PrincipalMax = prinstress1 - results.PrincipalMed = prinstress2 - results.PrincipalMin = prinstress3 - results.MaxShear = shearstress - - 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))) - results.NodeNumbers = stress.keys() - - # Read Equivalent Plastic strain if they exist - if 'peeq' in result_set: - Peeq = result_set['peeq'] - if len(Peeq) > 0: - if len(Peeq.values()) != len(disp.values()): - Pe = [] - Pe_extra_nodes = Peeq.values() - nodes = len(disp.values()) - for i in range(nodes): - Pe_value = Pe_extra_nodes[i] - Pe.append(Pe_value) - results.Peeq = Pe - else: - results.Peeq = Peeq.values() - - x_min, y_min, z_min = map(min, zip(*displacement)) - sum_list = map(sum, zip(*displacement)) - x_avg, y_avg, z_avg = [i / no_of_values for i in sum_list] - - s_max = max(results.StressValues) - s_min = min(results.StressValues) - s_avg = sum(results.StressValues) / no_of_values - p1_min = min(results.PrincipalMax) - p1_avg = sum(results.PrincipalMax) / no_of_values - p1_max = max(results.PrincipalMax) - p2_min = min(results.PrincipalMed) - p2_avg = sum(results.PrincipalMed) / no_of_values - p2_max = max(results.PrincipalMed) - p3_min = min(results.PrincipalMin) - p3_avg = sum(results.PrincipalMin) / no_of_values - p3_max = max(results.PrincipalMin) - ms_min = min(results.MaxShear) - ms_avg = sum(results.MaxShear) / no_of_values - ms_max = max(results.MaxShear) - if results.Peeq: - peeq_max = max(results.Peeq) - peeq_min = min(results.Peeq) - peeq_avg = sum(results.Peeq) / no_of_values - else: - peeq_max = 0 - peeq_min = 0 - peeq_avg = 0 - - 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) / no_of_values - - 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, - p1_min, p1_avg, p1_max, - p2_min, p2_avg, p2_max, - p3_min, p3_avg, p3_max, - ms_min, ms_avg, ms_max, - peeq_min, peeq_avg, peeq_max] - except: - pass - - # Read temperatures if they exist - try: - Temperature = result_set['temp'] - if len(Temperature) > 0: - if len(Temperature.values()) != len(disp.values()): - Temp = [] - Temp_extra_nodes = Temperature.values() - nodes = len(disp.values()) - for i in range(nodes): - Temp_value = Temp_extra_nodes[i] - Temp.append(Temp_value) - results.Temperature = list(map((lambda x: x), Temp)) - else: - results.Temperature = list(map((lambda x: x), Temperature.values())) - results.Time = step_time - except: - pass - - try: - MassFlow = result_set['mflow'] - if len(MassFlow) > 0: - results.MassFlowRate = list(map((lambda x: x), MassFlow.values())) - results.Time = step_time - except: - pass - - try: - NetworkPressure = result_set['npressure'] - if len(NetworkPressure) > 0: - results.NetworkPressure = list(map((lambda x: x), NetworkPressure.values())) - results.Time = step_time - except: - pass - + results = importToolsFem.fill_femresult_mechanical(results, result_set, span) analysis_object.Member = analysis_object.Member + [results] if(FreeCAD.GuiUp): @@ -690,32 +530,3 @@ def readResult(frd_input): 'Penta15Elem': elements_penta15, 'Hexa20Elem': elements_hexa20, 'Tria3Elem': elements_tria3, 'Tria6Elem': elements_tria6, 'Quad4Elem': elements_quad4, 'Quad8Elem': elements_quad8, 'Seg2Elem': elements_seg2, 'Seg3Elem': elements_seg3, 'Results': results} - - -# helper -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 calculate_principal_stress(i): - sigma = np.array([[i[0], i[3], i[4]], - [i[3], i[1], i[5]], - [i[4], i[5], i[2]]]) - # compute principal stresses - eigvals = list(np.linalg.eigvalsh(sigma)) - eigvals.sort() - eigvals.reverse() - maxshear = (eigvals[0] - eigvals[2]) / 2.0 - return (eigvals[0], eigvals[1], eigvals[2], maxshear) diff --git a/src/Mod/Fem/importToolsFem.py b/src/Mod/Fem/importToolsFem.py index 4ef3d655e..2fd91df41 100644 --- a/src/Mod/Fem/importToolsFem.py +++ b/src/Mod/Fem/importToolsFem.py @@ -29,6 +29,8 @@ __url__ = "http://www.freecadweb.org" # \brief FreeCAD FEM import tools import FreeCAD +from math import pow, sqrt +import numpy as np def make_femmesh(mesh_data): @@ -111,3 +113,200 @@ def make_femmesh(mesh_data): else: FreeCAD.Console.PrintError("No Nodes found!\n") return mesh + + +def fill_femresult_mechanical(results, result_set, span): + ''' fills an FreeCAD FEM mechanical result object with result data + ''' + + eigenmode_number = result_set['number'] + step_time = result_set['time'] + step_time = round(step_time, 2) + + try: + disp = result_set['disp'] + stressv = result_set['stressv'] + strainv = result_set['strainv'] + no_of_values = len(disp) + displacement = [] + for k, v in disp.items(): + displacement.append(v) + + x_max, y_max, z_max = map(max, zip(*displacement)) + if eigenmode_number > 0: + 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 = list(map((lambda x: x * scale), disp.values())) + results.StressVectors = list(map((lambda x: x * scale), stressv.values())) + results.StrainVectors = list(map((lambda x: x * scale), strainv.values())) + results.NodeNumbers = disp.keys() + + stress = result_set['stress'] + if len(stress) > 0: + mstress = [] + prinstress1 = [] + prinstress2 = [] + prinstress3 = [] + shearstress = [] + for i in stress.values(): + mstress.append(calculate_von_mises(i)) + prin1, prin2, prin3, shear = calculate_principal_stress(i) + prinstress1.append(prin1) + prinstress2.append(prin2) + prinstress3.append(prin3) + shearstress.append(shear) + if eigenmode_number > 0: + results.StressValues = list(map((lambda x: x * scale), mstress)) + results.PrincipalMax = list(map((lambda x: x * scale), prinstress1)) + results.PrincipalMed = list(map((lambda x: x * scale), prinstress2)) + results.PrincipalMin = list(map((lambda x: x * scale), prinstress3)) + results.MaxShear = list(map((lambda x: x * scale), shearstress)) + results.Eigenmode = eigenmode_number + else: + results.StressValues = mstress + results.PrincipalMax = prinstress1 + results.PrincipalMed = prinstress2 + results.PrincipalMin = prinstress3 + results.MaxShear = shearstress + + 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))) + results.NodeNumbers = stress.keys() + + # Read Equivalent Plastic strain if they exist + if 'peeq' in result_set: + Peeq = result_set['peeq'] + if len(Peeq) > 0: + if len(Peeq.values()) != len(disp.values()): + Pe = [] + Pe_extra_nodes = Peeq.values() + nodes = len(disp.values()) + for i in range(nodes): + Pe_value = Pe_extra_nodes[i] + Pe.append(Pe_value) + results.Peeq = Pe + else: + results.Peeq = Peeq.values() + + x_min, y_min, z_min = map(min, zip(*displacement)) + sum_list = map(sum, zip(*displacement)) + x_avg, y_avg, z_avg = [i / no_of_values for i in sum_list] + + s_max = max(results.StressValues) + s_min = min(results.StressValues) + s_avg = sum(results.StressValues) / no_of_values + p1_min = min(results.PrincipalMax) + p1_avg = sum(results.PrincipalMax) / no_of_values + p1_max = max(results.PrincipalMax) + p2_min = min(results.PrincipalMed) + p2_avg = sum(results.PrincipalMed) / no_of_values + p2_max = max(results.PrincipalMed) + p3_min = min(results.PrincipalMin) + p3_avg = sum(results.PrincipalMin) / no_of_values + p3_max = max(results.PrincipalMin) + ms_min = min(results.MaxShear) + ms_avg = sum(results.MaxShear) / no_of_values + ms_max = max(results.MaxShear) + if results.Peeq: + peeq_max = max(results.Peeq) + peeq_min = min(results.Peeq) + peeq_avg = sum(results.Peeq) / no_of_values + else: + peeq_max = 0 + peeq_min = 0 + peeq_avg = 0 + + 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) / no_of_values + + 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, + p1_min, p1_avg, p1_max, + p2_min, p2_avg, p2_max, + p3_min, p3_avg, p3_max, + ms_min, ms_avg, ms_max, + peeq_min, peeq_avg, peeq_max] + except: + pass + + # Read temperatures if they exist + try: + Temperature = result_set['temp'] + if len(Temperature) > 0: + if len(Temperature.values()) != len(disp.values()): + Temp = [] + Temp_extra_nodes = Temperature.values() + nodes = len(disp.values()) + for i in range(nodes): + Temp_value = Temp_extra_nodes[i] + Temp.append(Temp_value) + results.Temperature = list(map((lambda x: x), Temp)) + else: + results.Temperature = list(map((lambda x: x), Temperature.values())) + results.Time = step_time + except: + pass + + try: + MassFlow = result_set['mflow'] + if len(MassFlow) > 0: + results.MassFlowRate = list(map((lambda x: x), MassFlow.values())) + results.Time = step_time + except: + pass + + try: + NetworkPressure = result_set['npressure'] + if len(NetworkPressure) > 0: + results.NetworkPressure = list(map((lambda x: x), NetworkPressure.values())) + results.Time = step_time + except: + pass + + return results + + +# helper +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 calculate_principal_stress(i): + sigma = np.array([[i[0], i[3], i[4]], + [i[3], i[1], i[5]], + [i[4], i[5], i[2]]]) + # compute principal stresses + eigvals = list(np.linalg.eigvalsh(sigma)) + eigvals.sort() + eigvals.reverse() + maxshear = (eigvals[0] - eigvals[2]) / 2.0 + return (eigvals[0], eigvals[1], eigvals[2], maxshear)