diff --git a/src/Mod/Fem/feminout/importCcxFrdResults.py b/src/Mod/Fem/feminout/importCcxFrdResults.py index a6572232c2..3d9061d44c 100644 --- a/src/Mod/Fem/feminout/importCcxFrdResults.py +++ b/src/Mod/Fem/feminout/importCcxFrdResults.py @@ -105,6 +105,7 @@ def importFrd(filename, analysis=None, result_name_prefix=None): res_obj = restools.compact_result(res_obj) res_obj = restools.add_disp_apps(res_obj) # fill DisplacementLengths res_obj = restools.add_von_mises(res_obj) # fill StressValues + res_obj = restools.add_principal_stress(res_obj) # fill PrincipalMax, PrincipalMed, PrincipalMin, MaxShear res_obj = restools.fill_femresult_stats(res_obj) # fill Stats else: error_message = ( diff --git a/src/Mod/Fem/feminout/importToolsFem.py b/src/Mod/Fem/feminout/importToolsFem.py index 6bb6591bbd..71fc0a09ae 100644 --- a/src/Mod/Fem/feminout/importToolsFem.py +++ b/src/Mod/Fem/feminout/importToolsFem.py @@ -29,8 +29,6 @@ __url__ = "http://www.freecadweb.org" # \brief FreeCAD FEM import tools import FreeCAD -from math import pow, sqrt -import numpy as np def get_FemMeshObjectMeshGroups(fem_mesh_obj): @@ -310,32 +308,6 @@ def fill_femresult_mechanical(res_obj, result_set): strainv[i] = (FreeCAD.Vector(values_E[0], values_E[1], values_E[2])) res_obj.StrainVectors = list(map((lambda x: x * scale), strainv.values())) - # calculate von Mises, principal and max Shear and fill them in res_obj - if 'stress' in result_set: - stress = result_set['stress'] - if len(stress) > 0: - prinstress1 = [] - prinstress2 = [] - prinstress3 = [] - shearstress = [] - for i in stress.values(): - prin1, prin2, prin3, shear = calculate_principal_stress(i) - prinstress1.append(prin1) - prinstress2.append(prin2) - prinstress3.append(prin3) - shearstress.append(shear) - if eigenmode_number > 0: - res_obj.PrincipalMax = list(map((lambda x: x * scale), prinstress1)) - res_obj.PrincipalMed = list(map((lambda x: x * scale), prinstress2)) - res_obj.PrincipalMin = list(map((lambda x: x * scale), prinstress3)) - res_obj.MaxShear = list(map((lambda x: x * scale), shearstress)) - res_obj.Eigenmode = eigenmode_number - else: - res_obj.PrincipalMax = prinstress1 - res_obj.PrincipalMed = prinstress2 - res_obj.PrincipalMin = prinstress3 - res_obj.MaxShear = shearstress - # fill Equivalent Plastic strain if they exist if 'peeq' in result_set: Peeq = result_set['peeq'] @@ -387,23 +359,6 @@ def fill_femresult_mechanical(res_obj, result_set): # helper -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]]]) # https://forum.freecadweb.org/viewtopic.php?f=18&t=24637&start=10#p240408 - - try: # it will fail if NaN is inside the array, - # 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) - except: - return (float('NaN'), float('NaN'), float('NaN'), float('NaN')) - # TODO might be possible without a try except for NaN, https://forum.freecadweb.org/viewtopic.php?f=22&t=33911&start=10#p284229 - - def get_span(node_items): positions = [] # list of node vectors for k, v in node_items: diff --git a/src/Mod/Fem/femresult/resulttools.py b/src/Mod/Fem/femresult/resulttools.py index 784ed27a7e..e8d24e328e 100644 --- a/src/Mod/Fem/femresult/resulttools.py +++ b/src/Mod/Fem/femresult/resulttools.py @@ -279,6 +279,33 @@ def add_von_mises(res_obj): return res_obj +def add_principal_stress(res_obj): + prinstress1 = [] + prinstress2 = [] + prinstress3 = [] + shearstress = [] + iterator = zip( + res_obj.NodeStressXX, + res_obj.NodeStressYY, + res_obj.NodeStressZZ, + res_obj.NodeStressXY, + res_obj.NodeStressXZ, + res_obj.NodeStressYZ + ) + for Sxx, Syy, Szz, Sxy, Sxz, Syz in iterator: + prin1, prin2, prin3, shear = calculate_principal_stress((Sxx, Syy, Szz, Sxy, Sxz, Syz)) + prinstress1.append(prin1) + prinstress2.append(prin2) + prinstress3.append(prin3) + shearstress.append(shear) + res_obj.PrincipalMax = prinstress1 + res_obj.PrincipalMed = prinstress2 + res_obj.PrincipalMin = prinstress3 + res_obj.MaxShear = shearstress + FreeCAD.Console.PrintMessage('Added principal stress and max shear values.\n') + return res_obj + + def compact_result(res_obj): ''' compacts result.Mesh and appropriate result.NodeNumbers @@ -320,6 +347,24 @@ def calculate_von_mises(i): return vm_stress +def calculate_principal_stress(i): + import numpy as np + sigma = np.array([[i[0], i[3], i[4]], + [i[3], i[1], i[5]], + [i[4], i[5], i[2]]]) # https://forum.freecadweb.org/viewtopic.php?f=18&t=24637&start=10#p240408 + + try: # it will fail if NaN is inside the array, + # 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) + except: + return (float('NaN'), float('NaN'), float('NaN'), float('NaN')) + # TODO might be possible without a try except for NaN, https://forum.freecadweb.org/viewtopic.php?f=22&t=33911&start=10#p284229 + + def calculate_disp_abs(displacements): disp_abs = [] for d in displacements: