FEM: import tools, move fill result obj with imported results inside import tools module
This commit is contained in:
Bernd Hahnebach
wmayer
parent
a569319510
commit
6f105c4b67
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user