439 lines
18 KiB
Python
439 lines
18 KiB
Python
# ***************************************************************************
|
|
# * *
|
|
# * Copyright (c) 2017 - Bernd Hahnebach <bernd@bimstatik.org> *
|
|
# * *
|
|
# * This program is free software; you can redistribute it and/or modify *
|
|
# * it under the terms of the GNU Lesser General Public License (LGPL) *
|
|
# * as published by the Free Software Foundation; either version 2 of *
|
|
# * the License, or (at your option) any later version. *
|
|
# * for detail see the LICENCE text file. *
|
|
# * *
|
|
# * This program is distributed in the hope that it will be useful, *
|
|
# * but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
# * GNU Library General Public License for more details. *
|
|
# * *
|
|
# * You should have received a copy of the GNU Library General Public *
|
|
# * License along with this program; if not, write to the Free Software *
|
|
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 *
|
|
# * USA *
|
|
# * *
|
|
# ***************************************************************************
|
|
|
|
__title__ = "FreeCAD FEM import tools"
|
|
__author__ = "Bernd Hahnebach"
|
|
__url__ = "http://www.freecadweb.org"
|
|
|
|
## @package importToolsFem
|
|
# \ingroup FEM
|
|
# \brief FreeCAD FEM import tools
|
|
|
|
import FreeCAD
|
|
from math import pow, sqrt
|
|
import numpy as np
|
|
|
|
|
|
def get_FemMeshObjectMeshGroups(fem_mesh_obj):
|
|
"""
|
|
Get mesh groups from mesh. This also throws no exception if there
|
|
is no Groups property at all (e.g. Netgen meshes).
|
|
"""
|
|
fem_mesh = fem_mesh_obj.FemMesh
|
|
try:
|
|
gmshgroups = fem_mesh.Groups
|
|
except:
|
|
gmshgroups = ()
|
|
|
|
return gmshgroups
|
|
|
|
|
|
def get_FemMeshObjectOrder(fem_mesh_obj):
|
|
"""
|
|
Gets element order. Element order counting based on number of nodes on
|
|
edges. Edge with 2 nodes -> linear elements, Edge with 3 nodes ->
|
|
quadratic elements, and so on. No edges in mesh -> not determined.
|
|
(Is this possible? Seems to be a very degenerate case.)
|
|
If there are edges with different number of nodes appearing, return
|
|
list of orders.
|
|
"""
|
|
presumable_order = None
|
|
|
|
edges = fem_mesh_obj.FemMesh.Edges
|
|
|
|
if edges != ():
|
|
edges_length_set = list({len(fem_mesh_obj.FemMesh.getElementNodes(e)) for e in edges})
|
|
# only need set to eliminate double entries
|
|
|
|
if len(edges_length_set) == 1:
|
|
presumable_order = edges_length_set[0] - 1
|
|
else:
|
|
presumable_order = [el - 1 for el in edges_length_set]
|
|
else:
|
|
print("Found no edges in mesh: Element order determination does not work without them.")
|
|
|
|
return presumable_order
|
|
|
|
|
|
def get_FemMeshObjectDimension(fem_mesh_obj):
|
|
""" Count all entities in an abstract sense, to distinguish which dimension the mesh is
|
|
(i.e. linemesh, facemesh, volumemesh)
|
|
"""
|
|
dim = None
|
|
|
|
if fem_mesh_obj.FemMesh.Nodes != ():
|
|
dim = 0
|
|
if fem_mesh_obj.FemMesh.Edges != ():
|
|
dim = 1
|
|
if fem_mesh_obj.FemMesh.Faces != ():
|
|
dim = 2
|
|
if fem_mesh_obj.FemMesh.Volumes != ():
|
|
dim = 3
|
|
|
|
return dim
|
|
|
|
|
|
def get_FemMeshObjectElementTypes(fem_mesh_obj, remove_zero_element_entries=True):
|
|
"""
|
|
Spit out all elements in the mesh with their appropriate dimension.
|
|
"""
|
|
FreeCAD_element_names_dims = {
|
|
"Node": 0, "Edge": 1, "Hexa": 3, "Polygon": 2, "Polyhedron": 3,
|
|
"Prism": 3, "Pyramid": 3, "Quadrangle": 2, "Tetra": 3, "Triangle": 2}
|
|
|
|
elements_list_with_zero = [(eval("fem_mesh_obj.FemMesh." + s + "Count"), s, d) for (s, d) in FreeCAD_element_names_dims.iteritems()]
|
|
# ugly but necessary
|
|
if remove_zero_element_entries:
|
|
elements_list = [(num, s, d) for (num, s, d) in elements_list_with_zero if num > 0]
|
|
else:
|
|
elements_list = elements_list_with_zero
|
|
|
|
return elements_list
|
|
|
|
|
|
def get_MaxDimElementFromList(elem_list):
|
|
"""
|
|
Gets element with the maximal dimension in the mesh to determine cells.
|
|
"""
|
|
elem_list.sort(key=lambda t: t[2])
|
|
return elem_list[-1]
|
|
|
|
|
|
def make_femmesh(mesh_data):
|
|
''' makes an FreeCAD FEM Mesh object from FEM Mesh data
|
|
'''
|
|
import Fem
|
|
mesh = Fem.FemMesh()
|
|
m = mesh_data
|
|
if ('Nodes' in m) and (len(m['Nodes']) > 0):
|
|
print("Found: nodes")
|
|
if (
|
|
('Seg2Elem' in m) or
|
|
('Seg3Elem' in m) or
|
|
('Tria3Elem' in m) or
|
|
('Tria6Elem' in m) or
|
|
('Quad4Elem' in m) or
|
|
('Quad8Elem' in m) or
|
|
('Tetra4Elem' in m) or
|
|
('Tetra10Elem' in m) or
|
|
('Penta6Elem' in m) or
|
|
('Penta15Elem' in m) or
|
|
('Hexa8Elem' in m) or
|
|
('Hexa20Elem' in m)
|
|
):
|
|
|
|
nds = m['Nodes']
|
|
print("Found: elements")
|
|
for i in nds:
|
|
n = nds[i]
|
|
mesh.addNode(n[0], n[1], n[2], i)
|
|
elms_hexa8 = m['Hexa8Elem']
|
|
for i in elms_hexa8:
|
|
e = elms_hexa8[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
|
|
elms_penta6 = m['Penta6Elem']
|
|
for i in elms_penta6:
|
|
e = elms_penta6[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5]], i)
|
|
elms_tetra4 = m['Tetra4Elem']
|
|
for i in elms_tetra4:
|
|
e = elms_tetra4[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3]], i)
|
|
elms_tetra10 = m['Tetra10Elem']
|
|
for i in elms_tetra10:
|
|
e = elms_tetra10[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9]], i)
|
|
elms_penta15 = m['Penta15Elem']
|
|
for i in elms_penta15:
|
|
e = elms_penta15[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
|
|
e[10], e[11], e[12], e[13], e[14]], i)
|
|
elms_hexa20 = m['Hexa20Elem']
|
|
for i in elms_hexa20:
|
|
e = elms_hexa20[i]
|
|
mesh.addVolume([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7], e[8], e[9],
|
|
e[10], e[11], e[12], e[13], e[14], e[15], e[16], e[17], e[18], e[19]], i)
|
|
elms_tria3 = m['Tria3Elem']
|
|
for i in elms_tria3:
|
|
e = elms_tria3[i]
|
|
mesh.addFace([e[0], e[1], e[2]], i)
|
|
elms_tria6 = m['Tria6Elem']
|
|
for i in elms_tria6:
|
|
e = elms_tria6[i]
|
|
mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5]], i)
|
|
elms_quad4 = m['Quad4Elem']
|
|
for i in elms_quad4:
|
|
e = elms_quad4[i]
|
|
mesh.addFace([e[0], e[1], e[2], e[3]], i)
|
|
elms_quad8 = m['Quad8Elem']
|
|
for i in elms_quad8:
|
|
e = elms_quad8[i]
|
|
mesh.addFace([e[0], e[1], e[2], e[3], e[4], e[5], e[6], e[7]], i)
|
|
elms_seg2 = m['Seg2Elem']
|
|
for i in elms_seg2:
|
|
e = elms_seg2[i]
|
|
mesh.addEdge([e[0], e[1]], i)
|
|
elms_seg3 = m['Seg3Elem']
|
|
for i in elms_seg3:
|
|
e = elms_seg3[i]
|
|
mesh.addEdge([e[0], e[1], e[2]], i)
|
|
print("imported mesh: {} nodes, {} HEXA8, {} PENTA6, {} TETRA4, {} TETRA10, {} PENTA15".format(
|
|
len(nds), len(elms_hexa8), len(elms_penta6), len(elms_tetra4), len(elms_tetra10), len(elms_penta15)))
|
|
print("imported mesh: {} HEXA20, {} TRIA3, {} TRIA6, {} QUAD4, {} QUAD8, {} SEG2, {} SEG3".format(
|
|
len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2), len(elms_seg3)))
|
|
else:
|
|
FreeCAD.Console.PrintError("No Elements found!\n")
|
|
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
|
|
'''
|
|
no_of_values = None
|
|
|
|
if 'number' in result_set:
|
|
eigenmode_number = result_set['number']
|
|
else:
|
|
eigenmode_number = 0
|
|
if 'time' in result_set:
|
|
step_time = result_set['time']
|
|
step_time = round(step_time, 2)
|
|
|
|
if 'disp' in result_set:
|
|
disp = result_set['disp']
|
|
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
|
|
|
|
results.DisplacementVectors = list(map((lambda x: x * scale), disp.values()))
|
|
results.NodeNumbers = list(disp.keys())
|
|
results.DisplacementLengths = calculate_disp_abs(displacement)
|
|
|
|
if 'stressv' in result_set:
|
|
stressv = result_set['stressv']
|
|
results.StressVectors = list(map((lambda x: x * scale), stressv.values()))
|
|
|
|
if 'strainv' in result_set:
|
|
strainv = result_set['strainv']
|
|
results.StrainVectors = list(map((lambda x: x * scale), strainv.values()))
|
|
|
|
if 'stress' in result_set:
|
|
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
|
|
stress_keys = list(stress.keys())
|
|
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()
|
|
|
|
# Read temperatures if they exist
|
|
if 'temp' in result_set:
|
|
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
|
|
|
|
# read MassFlow, disp does not exist, no_of_values and results.NodeNumbers needs to be set
|
|
if 'mflow' in result_set:
|
|
MassFlow = result_set['mflow']
|
|
if len(MassFlow) > 0:
|
|
results.MassFlowRate = list(map((lambda x: x), MassFlow.values()))
|
|
results.Time = step_time
|
|
no_of_values = len(MassFlow)
|
|
results.NodeNumbers = list(MassFlow.keys())
|
|
|
|
# read NetworkPressure, disp does not exist, see MassFlow
|
|
if 'npressure' in result_set:
|
|
NetworkPressure = result_set['npressure']
|
|
if len(NetworkPressure) > 0:
|
|
results.NetworkPressure = list(map((lambda x: x), NetworkPressure.values()))
|
|
results.Time = step_time
|
|
|
|
# result stats, set stats values to 0, they may not exist
|
|
x_min = y_min = z_min = x_max = y_max = z_max = x_avg = y_avg = z_avg = 0
|
|
a_max = a_min = a_avg = s_max = s_min = s_avg = 0
|
|
p1_min = p1_avg = p1_max = p2_min = p2_avg = p2_max = p3_min = p3_avg = p3_max = 0
|
|
ms_min = ms_avg = ms_max = peeq_min = peeq_avg = peeq_max = 0
|
|
temp_min = temp_avg = temp_max = mflow_min = mflow_avg = mflow_max = npress_min = npress_avg = npress_max = 0
|
|
|
|
if results.DisplacementVectors:
|
|
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 / no_of_values for i in sum_list]
|
|
a_min = min(results.DisplacementLengths)
|
|
a_avg = sum(results.DisplacementLengths) / no_of_values
|
|
a_max = max(results.DisplacementLengths)
|
|
if results.StressValues:
|
|
s_min = min(results.StressValues)
|
|
s_avg = sum(results.StressValues) / no_of_values
|
|
s_max = max(results.StressValues)
|
|
if results.PrincipalMax:
|
|
p1_min = min(results.PrincipalMax)
|
|
p1_avg = sum(results.PrincipalMax) / no_of_values
|
|
p1_max = max(results.PrincipalMax)
|
|
if results.PrincipalMed:
|
|
p2_min = min(results.PrincipalMed)
|
|
p2_avg = sum(results.PrincipalMed) / no_of_values
|
|
p2_max = max(results.PrincipalMed)
|
|
if results.PrincipalMin:
|
|
p3_min = min(results.PrincipalMin)
|
|
p3_avg = sum(results.PrincipalMin) / no_of_values
|
|
p3_max = max(results.PrincipalMin)
|
|
if results.MaxShear:
|
|
ms_min = min(results.MaxShear)
|
|
ms_avg = sum(results.MaxShear) / no_of_values
|
|
ms_max = max(results.MaxShear)
|
|
if results.Peeq:
|
|
peeq_min = min(results.Peeq)
|
|
peeq_avg = sum(results.Peeq) / no_of_values
|
|
peeq_max = max(results.Peeq)
|
|
if results.Temperature:
|
|
temp_min = min(results.Temperature)
|
|
temp_avg = sum(results.Temperature) / no_of_values
|
|
temp_max = max(results.Temperature)
|
|
if results.MassFlowRate:
|
|
mflow_min = min(results.MassFlowRate)
|
|
mflow_avg = sum(results.MassFlowRate) / no_of_values
|
|
mflow_max = max(results.MassFlowRate)
|
|
if results.NetworkPressure:
|
|
npress_min = min(results.NetworkPressure)
|
|
npress_avg = sum(results.NetworkPressure) / no_of_values
|
|
npress_max = max(results.NetworkPressure)
|
|
|
|
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,
|
|
temp_min, temp_avg, temp_max,
|
|
mflow_min, mflow_avg, mflow_max,
|
|
npress_min, npress_avg, npress_max]
|
|
# do not forget to adapt the def get_stats in FemResultTools module as well as the TestFem module
|
|
# stat_types = ["U1", "U2", "U3", "Uabs", "Sabs", "MaxPrin", "MidPrin", "MinPrin", "MaxShear", "Peeq", "Temp", "MFlow", "NPress"]
|
|
# TODO: a dictionary would be far robust than a list, but needs adapten in VTK too because of VTK result import
|
|
|
|
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)
|
|
|
|
|
|
def calculate_disp_abs(displacements):
|
|
disp_abs = []
|
|
for d in displacements:
|
|
disp_abs.append(sqrt(pow(d[0], 2) + pow(d[1], 2) + pow(d[2], 2)))
|
|
return disp_abs
|