1489 lines
68 KiB
Python
1489 lines
68 KiB
Python
# ***************************************************************************
|
|
# * *
|
|
# * Copyright (c) 2016 - 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__ = "Tools for the work with FEM meshes"
|
|
__author__ = "Bernd Hahnebach"
|
|
__url__ = "http://www.freecadweb.org"
|
|
|
|
## \addtogroup FEM
|
|
# @{
|
|
|
|
import FreeCAD
|
|
|
|
|
|
def get_femnodes_by_femobj_with_references(femmesh, femobj):
|
|
node_set = []
|
|
if femmesh.GroupCount:
|
|
node_set = get_femnode_set_from_group_data(femmesh, femobj)
|
|
# print('node_set_group: ', node_set)
|
|
if not node_set:
|
|
node_set = get_femnodes_by_references(femmesh, femobj['Object'].References)
|
|
# print('node_set_nogroup: ', node_set)
|
|
return node_set
|
|
|
|
|
|
def get_femelements_by_references(femmesh, femelement_table, references, femnodes_ele_table=None):
|
|
'''get the femelements for a list of references
|
|
'''
|
|
references_femelements = []
|
|
for ref in references:
|
|
ref_femnodes = get_femnodes_by_refshape(femmesh, ref) # femnodes for the current ref
|
|
if femnodes_ele_table:
|
|
# blind fast binary search, works for volumes only
|
|
references_femelements += get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, ref_femnodes) # femelements for all references
|
|
else:
|
|
# standars search
|
|
references_femelements += get_femelements_by_femnodes_std(femelement_table, ref_femnodes) # femelements for all references
|
|
return references_femelements
|
|
|
|
|
|
def get_femnodes_by_references(femmesh, references):
|
|
'''get the femnodes for a list of references
|
|
'''
|
|
references_femnodes = []
|
|
for ref in references:
|
|
references_femnodes += get_femnodes_by_refshape(femmesh, ref)
|
|
|
|
# return references_femnodes # keeps duplicate nodes, keeps node order
|
|
|
|
# if nodes are used for nodesets, duplicats should be removed
|
|
return list(set(references_femnodes)) # removes duplicate nodes, sortes node order
|
|
|
|
|
|
def get_femnodes_by_refshape(femmesh, ref):
|
|
nodes = []
|
|
for refelement in ref[1]:
|
|
r = get_element(ref[0], refelement) # the method getElement(element) does not return Solid elements
|
|
print(' ReferenceShape : ', r.ShapeType, ', ', ref[0].Name, ', ', ref[0].Label, ' --> ', refelement)
|
|
if r.ShapeType == 'Vertex':
|
|
nodes += femmesh.getNodesByVertex(r)
|
|
elif r.ShapeType == 'Edge':
|
|
nodes += femmesh.getNodesByEdge(r)
|
|
elif r.ShapeType == 'Face':
|
|
nodes += femmesh.getNodesByFace(r)
|
|
elif r.ShapeType == 'Solid':
|
|
nodes += femmesh.getNodesBySolid(r)
|
|
else:
|
|
print(' No Vertice, Edge, Face or Solid as reference shapes!')
|
|
return nodes
|
|
|
|
|
|
def get_femelement_table(femmesh):
|
|
""" get_femelement_table(femmesh): { elementid : [ nodeid, nodeid, ... , nodeid ] }"""
|
|
femelement_table = {}
|
|
if is_solid_femmesh(femmesh):
|
|
for i in femmesh.Volumes:
|
|
femelement_table[i] = femmesh.getElementNodes(i)
|
|
elif is_face_femmesh(femmesh):
|
|
for i in femmesh.Faces:
|
|
femelement_table[i] = femmesh.getElementNodes(i)
|
|
elif is_edge_femmesh(femmesh):
|
|
for i in femmesh.Edges:
|
|
femelement_table[i] = femmesh.getElementNodes(i)
|
|
else:
|
|
FreeCAD.Console.PrintError('Neither solid nor face nor edge femmesh!\n')
|
|
return femelement_table
|
|
|
|
|
|
def get_femnodes_ele_table(femnodes_mesh, femelement_table):
|
|
'''the femnodes_ele_table contains for each node its membership in elements
|
|
{nodeID : [[eleID, NodePosition], [], ...], nodeID : [[], [], ...], ...}
|
|
stored informatation are:
|
|
element number, the number of nodes per element, the position of the node in the element.
|
|
The position of the node in the element is coded as a set bit at that position in a bit array (integer)
|
|
Fixme: the number of nodes per element should be replaced by the type of the element
|
|
but I did not know, how to get this from the mesh.
|
|
Since the femelement_table contains either volume or face or edgemesh the femnodes_ele_table only
|
|
has either volume or face or edge elements, see get_femelement_table()
|
|
'''
|
|
femnodes_ele_table = {} # node_dict in ulrichs class
|
|
for n in femnodes_mesh: # initialize it with sorted node keys and empty lists
|
|
femnodes_ele_table[n] = []
|
|
for ele in femelement_table:
|
|
ele_list = femelement_table[ele]
|
|
# print(ele_list)
|
|
pos = int(1)
|
|
for ele_node in ele_list:
|
|
femnodes_ele_table[ele_node].append([ele, pos])
|
|
pos = pos << 1
|
|
print('len femnodes_ele_table:' + str(len(femnodes_ele_table)))
|
|
# print('femnodes_ele_table: ', femnodes_ele_table)
|
|
return femnodes_ele_table
|
|
|
|
|
|
def get_copy_of_empty_femelement_table(femelement_table):
|
|
'''{eleID : 0, eleID : 0, ...}
|
|
'''
|
|
empty_femelement_table = {}
|
|
for ele in femelement_table: # initialize it with sorted element keys and empty int
|
|
empty_femelement_table[ele] = 0
|
|
return empty_femelement_table.copy()
|
|
|
|
|
|
def get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_set):
|
|
'''Now we are looking for nodes inside of the Faces = filling the bit_pattern_dict
|
|
{eleID : [lenEleNodes, binary_position]}
|
|
see forumpost for a ver good explanation whats really happening
|
|
http://forum.freecadweb.org/viewtopic.php?f=18&p=141133&sid=013c93f496a63872951d2ce521702ffa#p141108
|
|
The bit_pattern_dict holds later an integer (bit array) for each element, which gives us
|
|
the information we are searching for:
|
|
Is this element part of the node list (searching for elements) or has this element a face we are searching for?
|
|
The number in the ele_dict is organized as a bit array.
|
|
The corresponding bit is set, if the node of the node_set is contained in the element.
|
|
'''
|
|
print('len femnodes_ele_table:' + str(len(femnodes_ele_table)))
|
|
print('len node_set: ' + str(len(node_set)))
|
|
# print('node_set: ', node_set)
|
|
bit_pattern_dict = get_copy_of_empty_femelement_table(femelement_table)
|
|
# # initializing the bit_pattern_dict
|
|
for ele in femelement_table:
|
|
len_ele = len(femelement_table[ele])
|
|
bit_pattern_dict[ele] = [len_ele, 0]
|
|
for node in node_set:
|
|
for nList in femnodes_ele_table[node]:
|
|
bit_pattern_dict[nList[0]][1] += nList[1]
|
|
print('len bit_pattern_dict:' + str(len(bit_pattern_dict)))
|
|
# print('bit_pattern_dict: ', bit_pattern_dict)
|
|
return bit_pattern_dict
|
|
|
|
|
|
def get_ccxelement_faces_from_binary_search(bit_pattern_dict):
|
|
'''get the CalculiX element face numbers
|
|
'''
|
|
tet10_mask = {
|
|
119: 1,
|
|
411: 2,
|
|
717: 3,
|
|
814: 4}
|
|
tet4_mask = {
|
|
7: 1,
|
|
11: 2,
|
|
13: 3,
|
|
14: 4}
|
|
hex8_mask = {
|
|
240: 1,
|
|
15: 2,
|
|
102: 3,
|
|
204: 4,
|
|
153: 5,
|
|
51: 6}
|
|
hex20_mask = {
|
|
61680: 1,
|
|
3855: 2,
|
|
402022: 3,
|
|
804044: 4,
|
|
624793: 5,
|
|
201011: 6}
|
|
pent6_mask = {
|
|
56: 1,
|
|
7: 2,
|
|
54: 3,
|
|
45: 4,
|
|
27: 5}
|
|
pent15_mask = {
|
|
3640: 1,
|
|
455: 2,
|
|
25782: 3,
|
|
22829: 4,
|
|
12891: 5}
|
|
vol_dict = {
|
|
4: tet4_mask,
|
|
6: pent6_mask,
|
|
8: hex8_mask,
|
|
10: tet10_mask,
|
|
15: pent15_mask,
|
|
20: hex20_mask}
|
|
faces = []
|
|
for ele in bit_pattern_dict:
|
|
mask_dict = vol_dict[bit_pattern_dict[ele][0]]
|
|
for key in mask_dict:
|
|
if (key & bit_pattern_dict[ele][1]) == key:
|
|
faces.append([ele, mask_dict[key]])
|
|
print('found Faces: ', len(faces))
|
|
print('faces: ', faces)
|
|
return faces
|
|
|
|
|
|
def get_femelements_by_femnodes_bin(femelement_table, femnodes_ele_table, node_list):
|
|
'''for every femelement of femelement_table
|
|
if all nodes of the femelement are in node_list,
|
|
the femelement is added to the list which is returned
|
|
blind fast binary search, but workd for volumes only
|
|
'''
|
|
print('binary search: get_femelements_by_femnodes_bin')
|
|
vol_masks = {
|
|
4: 15,
|
|
6: 63,
|
|
8: 255,
|
|
10: 1023,
|
|
15: 32767,
|
|
20: 1048575}
|
|
# Now we are looking for nodes inside of the Volumes = filling the bit_pattern_dict
|
|
print('len femnodes_ele_table:' + str(len(femnodes_ele_table)))
|
|
bit_pattern_dict = get_bit_pattern_dict(femelement_table, femnodes_ele_table, node_list)
|
|
# search
|
|
ele_list = [] # The ele_list contains the result of the search.
|
|
for ele in bit_pattern_dict:
|
|
# print('bit_pattern_dict[ele][0]: ', bit_pattern_dict[ele][0])
|
|
if bit_pattern_dict[ele][1] == vol_masks[bit_pattern_dict[ele][0]]:
|
|
ele_list.append(ele)
|
|
print('found Volumes: ', len(ele_list))
|
|
print(' volumes: ', len(ele_list))
|
|
return ele_list
|
|
|
|
|
|
def get_femelements_by_femnodes_std(femelement_table, node_list):
|
|
'''for every femelement of femelement_table
|
|
if all nodes of the femelement are in node_list,
|
|
the femelement is added to the list which is returned
|
|
e: elementlist
|
|
nodes: nodelist '''
|
|
print('std search: get_femelements_by_femnodes_std')
|
|
e = [] # elementlist
|
|
for elementID in sorted(femelement_table):
|
|
nodecount = 0
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == len(femelement_table[elementID]): # all nodes of the element are in the node_list!
|
|
e.append(elementID)
|
|
return e
|
|
|
|
|
|
def get_femvolumeelements_by_femfacenodes(femelement_table, node_list):
|
|
'''assume femelement_table only has volume elements
|
|
for every femvolumeelement of femelement_table
|
|
for tetra4 and tetra10 the C++ methods could be used --> test again to be sure
|
|
if hexa8 volume element --> if exact 4 element nodes are in node_list --> add femelement
|
|
if hexa20 volume element --> if exact 8 element nodes are in node_list --> add femelement
|
|
if penta6 volume element --> if exact 3 or 6 element nodes are in node_list --> add femelement
|
|
if penta15 volume element --> if exact 6 or 8 element nodes are in node_list --> add femelement
|
|
e: elementlist
|
|
nodes: nodelist '''
|
|
e = [] # elementlist
|
|
for elementID in sorted(femelement_table):
|
|
nodecount = 0
|
|
el_nd_ct = len(femelement_table[elementID])
|
|
if el_nd_ct == 4: # tetra4
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 3:
|
|
e.append(elementID)
|
|
elif el_nd_ct == 10: # tetra10
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 4:
|
|
e.append(elementID)
|
|
elif el_nd_ct == 8: # hexa8
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 4:
|
|
e.append(elementID)
|
|
elif el_nd_ct == 20: # hexa20
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 8:
|
|
e.append(elementID)
|
|
elif el_nd_ct == 6: # penta6
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 3 or nodecount == 4:
|
|
e.append(elementID)
|
|
elif el_nd_ct == 15: # penta15
|
|
for nodeID in femelement_table[elementID]:
|
|
if nodeID in node_list:
|
|
nodecount = nodecount + 1
|
|
if nodecount == 6 or nodecount == 8:
|
|
e.append(elementID)
|
|
else:
|
|
FreeCAD.Console.PrintError('Error in get_femvolumeelements_by_femfacenodes(): not known volume element: ' + el_nd_ct + '\n')
|
|
# print(sorted(e))
|
|
return e
|
|
|
|
|
|
def get_femelement_sets(femmesh, femelement_table, fem_objects, femnodes_ele_table=None): # fem_objects = FreeCAD FEM document objects
|
|
# get femelements for reference shapes of each obj.References
|
|
count_femelements = 0
|
|
referenced_femelements = []
|
|
has_remaining_femelements = None
|
|
for fem_object_i, fem_object in enumerate(fem_objects):
|
|
obj = fem_object['Object']
|
|
fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier
|
|
if obj.References:
|
|
ref_shape_femelements = []
|
|
ref_shape_femelements = get_femelements_by_references(femmesh, femelement_table, obj.References, femnodes_ele_table)
|
|
referenced_femelements += ref_shape_femelements
|
|
count_femelements += len(ref_shape_femelements)
|
|
fem_object['FEMElements'] = ref_shape_femelements
|
|
else:
|
|
has_remaining_femelements = obj.Name
|
|
# get remaining femelements for the fem_objects
|
|
if has_remaining_femelements:
|
|
remaining_femelements = []
|
|
for elemid in femelement_table:
|
|
if elemid not in referenced_femelements:
|
|
remaining_femelements.append(elemid)
|
|
count_femelements += len(remaining_femelements)
|
|
for fem_object in fem_objects:
|
|
obj = fem_object['Object']
|
|
if obj.Name == has_remaining_femelements:
|
|
fem_object['FEMElements'] = sorted(remaining_femelements)
|
|
# check if all worked out well
|
|
if not femelements_count_ok(len(femelement_table), count_femelements):
|
|
FreeCAD.Console.PrintError('Error in get_femelement_sets -- > femelements_count_ok() failed!\n')
|
|
|
|
|
|
def get_femnode_set_from_group_data(femmesh, fem_object):
|
|
# get femnodes from femmesh groupdata for reference shapes of each obj.References
|
|
# we assume the mesh group data fits with the reference shapes, no check is done in this regard !!!
|
|
# what happens if a reference shape was changed, but the mesh and the mesh groups were not created new !?!
|
|
obj = fem_object['Object']
|
|
group_nodes = None
|
|
if femmesh.GroupCount:
|
|
for g in femmesh.Groups:
|
|
grp_name = femmesh.getGroupName(g)
|
|
if grp_name.startswith(obj.Name):
|
|
if femmesh.getGroupElementType(g) == "Node":
|
|
print("Constraint: " + obj.Name + " --> " + "mesh group: " + grp_name)
|
|
group_nodes = femmesh.getGroupElements(g) # == ref_shape_femelements
|
|
return group_nodes
|
|
|
|
|
|
def get_femelement_sets_from_group_data(femmesh, fem_objects):
|
|
# get femelements from femmesh groupdata for reference shapes of each obj.References
|
|
count_femelements = 0
|
|
sum_group_elements = []
|
|
# we assume the mesh group data fits with the reference shapes, no check is done in this regard !!!
|
|
# what happens if a reference shape was changed, but the mesh and the mesh groups were not created new !?!
|
|
for fem_object_i, fem_object in enumerate(fem_objects):
|
|
obj = fem_object['Object']
|
|
fem_object['ShortName'] = get_elset_short_name(obj, fem_object_i) # unique short identifier
|
|
if femmesh.GroupCount:
|
|
for g in femmesh.Groups:
|
|
grp_name = femmesh.getGroupName(g)
|
|
if grp_name.startswith(obj.Name):
|
|
if femmesh.getGroupElementType(g) == "Volume":
|
|
print("Constraint: " + obj.Name + " --> " + "mesh group: " + grp_name)
|
|
group_elements = femmesh.getGroupElements(g) # == ref_shape_femelements
|
|
sum_group_elements += group_elements
|
|
count_femelements += len(group_elements)
|
|
fem_object['FEMElements'] = group_elements
|
|
# check if all worked out well
|
|
if not femelements_count_ok(femmesh.VolumeCount, count_femelements):
|
|
FreeCAD.Console.PrintError('Error in get_femelement_sets_from_group_data -- > femelements_count_ok() failed!\n')
|
|
return False
|
|
else:
|
|
return True
|
|
print("")
|
|
|
|
|
|
def get_elset_short_name(obj, i):
|
|
if hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemMaterial':
|
|
return 'Mat' + str(i)
|
|
elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemBeamSection':
|
|
return 'Beam' + str(i)
|
|
elif hasattr(obj, "Proxy") and obj.Proxy.Type == 'FemShellThickness':
|
|
return 'Shell' + str(i)
|
|
else:
|
|
print('Error: ', obj.Name, ' --> ', obj.Proxy.Type)
|
|
|
|
|
|
def get_force_obj_vertex_nodeload_table(femmesh, frc_obj):
|
|
# force_obj_node_load_table = [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)]
|
|
force_obj_node_load_table = []
|
|
node_load = frc_obj.Force / len(frc_obj.References)
|
|
for o, elem_tup in frc_obj.References:
|
|
node_count = len(elem_tup)
|
|
for elem in elem_tup:
|
|
ref_node = o.Shape.getElement(elem)
|
|
node = femmesh.getNodesByVertex(ref_node)
|
|
elem_info_string = 'node load on shape: ' + o.Name + ':' + elem
|
|
force_obj_node_load_table.append((elem_info_string, {node[0]: node_load / node_count}))
|
|
|
|
return force_obj_node_load_table
|
|
|
|
|
|
def get_force_obj_edge_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj):
|
|
# force_obj_node_load_table = [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)]
|
|
force_obj_node_load_table = []
|
|
sum_ref_edge_length = 0
|
|
sum_ref_edge_node_length = 0 # for debugging
|
|
sum_node_load = 0 # for debugging
|
|
for o, elem_tup in frc_obj.References:
|
|
for elem in elem_tup:
|
|
sum_ref_edge_length += o.Shape.getElement(elem).Length
|
|
if sum_ref_edge_length != 0:
|
|
force_per_sum_ref_edge_length = frc_obj.Force / sum_ref_edge_length
|
|
for o, elem_tup in frc_obj.References:
|
|
for elem in elem_tup:
|
|
ref_edge = o.Shape.getElement(elem)
|
|
|
|
# edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
|
|
edge_table = get_ref_edgenodes_table(femmesh, femelement_table, ref_edge)
|
|
|
|
# node_length_table = [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry
|
|
node_length_table = get_ref_edgenodes_lengths(femnodes_mesh, edge_table)
|
|
|
|
# node_sum_length_table = { nodeID : Length, ... , nodeID : Length } LengthSum for each node, one entry for each node
|
|
node_sum_length_table = get_ref_shape_node_sum_geom_table(node_length_table)
|
|
|
|
# node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad } NodeLoad for each node, one entry for each node
|
|
node_load_table = {}
|
|
sum_node_lengths = 0 # for debugging
|
|
for node in node_sum_length_table:
|
|
sum_node_lengths += node_sum_length_table[node] # for debugging
|
|
node_load_table[node] = node_sum_length_table[node] * force_per_sum_ref_edge_length
|
|
ratio_refedge_lengths = sum_node_lengths / ref_edge.Length
|
|
if ratio_refedge_lengths < 0.99 or ratio_refedge_lengths > 1.01:
|
|
FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
|
|
print(' sum_node_lengths:', sum_node_lengths)
|
|
print(' refedge_length: ', ref_edge.Length)
|
|
bad_refedge = ref_edge
|
|
sum_ref_edge_node_length += sum_node_lengths
|
|
|
|
elem_info_string = 'node loads on shape: ' + o.Name + ':' + elem
|
|
force_obj_node_load_table.append((elem_info_string, node_load_table))
|
|
|
|
for ref_shape in force_obj_node_load_table:
|
|
for node in ref_shape[1]:
|
|
sum_node_load += ref_shape[1][node] # for debugging
|
|
|
|
ratio = sum_node_load / frc_obj.Force
|
|
if ratio < 0.99 or ratio > 1.01:
|
|
print('Deviation sum_node_load to frc_obj.Force is more than 1% : ', ratio)
|
|
print(' sum_ref_edge_node_length: ', sum_ref_edge_node_length)
|
|
print(' sum_ref_edge_length: ', sum_ref_edge_length)
|
|
print(' sum_node_load: ', sum_node_load)
|
|
print(' frc_obj.Force: ', frc_obj.Force)
|
|
print(' the reason could be simply a circle length --> see method get_ref_edge_node_lengths')
|
|
print(' the reason could also be an problem in retrieving the ref_edge_node_length')
|
|
|
|
# try debugging of the last bad refedge
|
|
print('DEBUGGING')
|
|
print(bad_refedge)
|
|
|
|
print('bad_refedge_nodes')
|
|
bad_refedge_nodes = femmesh.getNodesByEdge(bad_refedge)
|
|
print(len(bad_refedge_nodes))
|
|
print(bad_refedge_nodes)
|
|
# import FreeCADGui
|
|
# FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_refedge_nodes
|
|
|
|
print('bad_edge_table')
|
|
# bad_edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
|
|
bad_edge_table = get_ref_edgenodes_table(femmesh, femelement_table, bad_refedge)
|
|
print(len(bad_edge_table))
|
|
bad_edge_table_nodes = []
|
|
for elem in bad_edge_table:
|
|
print(elem, ' --> ', bad_edge_table[elem])
|
|
for node in bad_edge_table[elem]:
|
|
if node not in bad_edge_table_nodes:
|
|
bad_edge_table_nodes.append(node)
|
|
print('sorted(bad_edge_table_nodes)')
|
|
print(sorted(bad_edge_table_nodes)) # should be == bad_refedge_nodes
|
|
# import FreeCADGui
|
|
# FreeCADGui.ActiveDocument.Compound_Mesh.HighlightedNodes = bad_edge_table_nodes
|
|
# bad_node_length_table = [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry
|
|
|
|
print('good_edge_table')
|
|
good_edge_table = delete_duplicate_mesh_elements(bad_edge_table)
|
|
for elem in good_edge_table:
|
|
print(elem, ' --> ', bad_edge_table[elem])
|
|
|
|
print('bad_node_length_table')
|
|
bad_node_length_table = get_ref_edgenodes_lengths(femnodes_mesh, bad_edge_table)
|
|
for n, l in bad_node_length_table:
|
|
print(n, ' --> ', l)
|
|
|
|
return force_obj_node_load_table
|
|
|
|
|
|
def get_pressure_obj_faces_depreciated(femmesh, femobj):
|
|
pressure_faces = []
|
|
for o, elem_tup in femobj['Object'].References:
|
|
for elem in elem_tup:
|
|
ref_shape = o.Shape.getElement(elem)
|
|
elem_info_string = 'face load on shape: ' + o.Name + ':' + elem
|
|
print(elem_info_string)
|
|
if ref_shape.ShapeType == 'Face':
|
|
pressure_faces.append((elem_info_string, femmesh.getccxVolumesByFace(ref_shape)))
|
|
return pressure_faces
|
|
|
|
|
|
def get_pressure_obj_faces(femmesh, femelement_table, femnodes_ele_table, femobj):
|
|
# get the nodes
|
|
prs_face_node_set = get_femnodes_by_femobj_with_references(femmesh, femobj) # sorted and duplicates removed
|
|
# print('prs_face_node_set: ', prs_face_node_set)
|
|
# fill the bit_pattern_dict and search for the faces
|
|
bit_pattern_dict = get_bit_pattern_dict(femelement_table, femnodes_ele_table, prs_face_node_set)
|
|
pressure_faces = get_ccxelement_faces_from_binary_search(bit_pattern_dict)
|
|
return pressure_faces
|
|
|
|
|
|
def get_force_obj_face_nodeload_table(femmesh, femelement_table, femnodes_mesh, frc_obj):
|
|
# force_obj_node_load_table = [('refshape_name.elemname',node_load_table), ..., ('refshape_name.elemname',node_load_table)]
|
|
force_obj_node_load_table = []
|
|
sum_ref_face_area = 0
|
|
sum_ref_face_node_area = 0 # for debugging
|
|
sum_node_load = 0 # for debugging
|
|
for o, elem_tup in frc_obj.References:
|
|
for elem in elem_tup:
|
|
sum_ref_face_area += o.Shape.getElement(elem).Area
|
|
if sum_ref_face_area != 0:
|
|
force_per_sum_ref_face_area = frc_obj.Force / sum_ref_face_area
|
|
for o, elem_tup in frc_obj.References:
|
|
for elem in elem_tup:
|
|
ref_face = o.Shape.getElement(elem)
|
|
|
|
# face_table = { meshfaceID : ( nodeID, ... , nodeID ) }
|
|
face_table = get_ref_facenodes_table(femmesh, femelement_table, ref_face)
|
|
|
|
# node_area_table = [ (nodeID, Area), ... , (nodeID, Area) ] some nodes will have more than one entry
|
|
node_area_table = get_ref_facenodes_areas(femnodes_mesh, face_table)
|
|
|
|
# node_sum_area_table = { nodeID : Area, ... , nodeID : Area } AreaSum for each node, one entry for each node
|
|
node_sum_area_table = get_ref_shape_node_sum_geom_table(node_area_table)
|
|
|
|
# node_load_table = { nodeID : NodeLoad, ... , nodeID : NodeLoad } NodeLoad for each node, one entry for each node
|
|
node_load_table = {}
|
|
sum_node_areas = 0 # for debugging
|
|
for node in node_sum_area_table:
|
|
sum_node_areas += node_sum_area_table[node] # for debugging
|
|
node_load_table[node] = node_sum_area_table[node] * force_per_sum_ref_face_area
|
|
ratio_refface_areas = sum_node_areas / ref_face.Area
|
|
if ratio_refface_areas < 0.99 or ratio_refface_areas > 1.01:
|
|
FreeCAD.Console.PrintError('Error on: ' + frc_obj.Name + ' --> ' + o.Name + '.' + elem + '\n')
|
|
print(' sum_node_areas:', sum_node_areas)
|
|
print(' ref_face_area: ', ref_face.Area)
|
|
sum_ref_face_node_area += sum_node_areas
|
|
|
|
elem_info_string = 'node loads on shape: ' + o.Name + ':' + elem
|
|
force_obj_node_load_table.append((elem_info_string, node_load_table))
|
|
|
|
for ref_shape in force_obj_node_load_table:
|
|
for node in ref_shape[1]:
|
|
sum_node_load += ref_shape[1][node] # for debugging
|
|
|
|
ratio = sum_node_load / frc_obj.Force
|
|
if ratio < 0.99 or ratio > 1.01:
|
|
print('Deviation sum_node_load to frc_obj.Force is more than 1% : ', ratio)
|
|
print(' sum_ref_face_node_area: ', sum_ref_face_node_area)
|
|
print(' sum_ref_face_area: ', sum_ref_face_area)
|
|
print(' sum_node_load: ', sum_node_load)
|
|
print(' frc_obj.Force: ', frc_obj.Force)
|
|
print(' the reason could be simply a circle area --> see method get_ref_face_node_areas')
|
|
print(' the reason could also be an problem in retrieving the ref_face_node_area')
|
|
|
|
return force_obj_node_load_table
|
|
|
|
|
|
def get_ref_edgenodes_table(femmesh, femelement_table, refedge):
|
|
edge_table = {} # { meshedgeID : ( nodeID, ... , nodeID ) }
|
|
refedge_nodes = femmesh.getNodesByEdge(refedge)
|
|
if is_solid_femmesh(femmesh):
|
|
refedge_fem_volumeelements = []
|
|
# if at least two nodes of a femvolumeelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
|
|
for elem in femelement_table:
|
|
nodecount = 0
|
|
for node in femelement_table[elem]:
|
|
if node in refedge_nodes:
|
|
nodecount += 1
|
|
if nodecount > 1:
|
|
refedge_fem_volumeelements.append(elem)
|
|
# for every refedge_fem_volumeelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
|
|
for elem in refedge_fem_volumeelements:
|
|
fe_refedge_nodes = []
|
|
for node in femelement_table[elem]:
|
|
if node in refedge_nodes:
|
|
fe_refedge_nodes.append(node)
|
|
edge_table[elem] = fe_refedge_nodes # { volumeID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
|
|
# FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
|
|
edge_table = delete_duplicate_mesh_elements(edge_table)
|
|
elif is_face_femmesh(femmesh):
|
|
refedge_fem_faceelements = []
|
|
# if at least two nodes of a femfaceelement are in refedge_nodes the volume is added to refedge_fem_volumeelements
|
|
for elem in femelement_table:
|
|
nodecount = 0
|
|
for node in femelement_table[elem]:
|
|
if node in refedge_nodes:
|
|
nodecount += 1
|
|
if nodecount > 1:
|
|
refedge_fem_faceelements.append(elem)
|
|
# for every refedge_fem_faceelement look which of his nodes is in refedge_nodes --> add all these nodes to edge_table
|
|
for elem in refedge_fem_faceelements:
|
|
fe_refedge_nodes = []
|
|
for node in femelement_table[elem]:
|
|
if node in refedge_nodes:
|
|
fe_refedge_nodes.append(node)
|
|
edge_table[elem] = fe_refedge_nodes # { faceID : ( edgenodeID, ... , edgenodeID )} # only the refedge nodes
|
|
# FIXME duplicate_mesh_elements: as soon as contact ans springs are supported the user should decide on which edge the load is applied
|
|
edge_table = delete_duplicate_mesh_elements(edge_table)
|
|
elif is_edge_femmesh(femmesh):
|
|
refedge_fem_edgeelements = get_femelements_by_femnodes_std(femelement_table, refedge_nodes)
|
|
for elem in refedge_fem_edgeelements:
|
|
edge_table[elem] = femelement_table[elem] # { edgeID : ( nodeID, ... , nodeID )} # all nodes off this femedgeelement
|
|
return edge_table
|
|
|
|
|
|
def get_ref_edgenodes_lengths(femnodes_mesh, edge_table):
|
|
# calulate the appropriate node_length for every node of every mesh edge (me)
|
|
# G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
|
|
|
|
# [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry
|
|
if (not femnodes_mesh) or (not edge_table):
|
|
FreeCAD.Console.PrintError("Error in get_ref_edgenodes_lengths(): Empty femnodes_mesh or edge_table!\n")
|
|
return []
|
|
node_length_table = []
|
|
mesh_edge_length = 0
|
|
# print(len(edge_table))
|
|
for me in edge_table:
|
|
femmesh_edgetype = len(edge_table[me])
|
|
if femmesh_edgetype == 2: # 2 node femmesh edge
|
|
# end_node_length = mesh_edge_length / 2
|
|
# ______
|
|
# P1 P2
|
|
P1 = femnodes_mesh[edge_table[me][0]]
|
|
P2 = femnodes_mesh[edge_table[me][1]]
|
|
edge_vec = P2 - P1
|
|
mesh_edge_length = edge_vec.Length
|
|
# print(mesh_edge_length)
|
|
end_node_length = mesh_edge_length / 2.0
|
|
node_length_table.append((edge_table[me][0], end_node_length))
|
|
node_length_table.append((edge_table[me][1], end_node_length))
|
|
|
|
elif femmesh_edgetype == 3: # 3 node femmesh edge
|
|
# end_node_length = mesh_edge_length / 6
|
|
# middle_node_length = mesh_face_area * 2 / 3
|
|
# _______ _______
|
|
# P1 P3 P2
|
|
P1 = femnodes_mesh[edge_table[me][0]]
|
|
P2 = femnodes_mesh[edge_table[me][1]]
|
|
P3 = femnodes_mesh[edge_table[me][2]]
|
|
edge_vec1 = P3 - P1
|
|
edge_vec2 = P2 - P3
|
|
mesh_edge_length = edge_vec1.Length + edge_vec2.Length
|
|
# print(me, ' --> ', mesh_edge_length)
|
|
end_node_length = mesh_edge_length / 6.0
|
|
middle_node_length = mesh_edge_length * 2.0 / 3.0
|
|
node_length_table.append((edge_table[me][0], end_node_length))
|
|
node_length_table.append((edge_table[me][1], end_node_length))
|
|
node_length_table.append((edge_table[me][2], middle_node_length))
|
|
return node_length_table
|
|
|
|
|
|
def get_ref_facenodes_table(femmesh, femelement_table, ref_face):
|
|
face_table = {} # { meshfaceID : ( nodeID, ... , nodeID ) }
|
|
if is_solid_femmesh(femmesh):
|
|
if has_no_face_data(femmesh):
|
|
print('No face date in volume mesh. We try to use getccxVolumesByFace() to retrive the volume elments of the ref_face!')
|
|
# there is no face data
|
|
# the problem if we retrive the nodes ourself is they are not sorted we just have the nodes. We need to sourt them according
|
|
# the shell mesh notaion of tria3, tria6, quad4, quad8
|
|
ref_face_nodes = femmesh.getNodesByFace(ref_face)
|
|
# try to use getccxVolumesByFace() to get the volume ids of element with elementfaces on the ref_face --> should work for tetra4 and tetra10
|
|
ref_face_volume_elements = femmesh.getccxVolumesByFace(ref_face) # list of tupels (mv, ccx_face_nr)
|
|
if ref_face_volume_elements: # mesh with tetras
|
|
print('Use of getccxVolumesByFace() has returned volume elements of the ref_face!')
|
|
for ve in ref_face_volume_elements:
|
|
veID = ve[0]
|
|
ve_ref_face_nodes = []
|
|
for nodeID in femelement_table[veID]:
|
|
if nodeID in ref_face_nodes:
|
|
ve_ref_face_nodes.append(nodeID)
|
|
face_table[veID] = ve_ref_face_nodes # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
|
|
else: # mesh with hexa or penta
|
|
print('Use of getccxVolumesByFace() has NOT returned volume elements of the ref_face! We try to use get_femvolumeelements_by_femfacenodes()!')
|
|
ref_face_volume_elements = get_femvolumeelements_by_femfacenodes(femelement_table, ref_face_nodes) # list of integer [mv]
|
|
for veID in ref_face_volume_elements:
|
|
ve_ref_face_nodes = []
|
|
for nodeID in femelement_table[veID]:
|
|
if nodeID in ref_face_nodes:
|
|
ve_ref_face_nodes.append(nodeID)
|
|
face_table[veID] = ve_ref_face_nodes # { volumeID : ( facenodeID, ... , facenodeID ) } only the ref_face nodes
|
|
face_table = build_mesh_faces_of_volume_elements(face_table, femelement_table) # we need to resort the nodes to make them build a element face
|
|
else: # the femmesh has face_data
|
|
faces = femmesh.getFacesByFace(ref_face) # (mv, mf)
|
|
for mf in faces:
|
|
face_table[mf] = femmesh.getElementNodes(mf)
|
|
elif is_face_femmesh(femmesh):
|
|
ref_face_nodes = femmesh.getNodesByFace(ref_face)
|
|
ref_face_elements = get_femelements_by_femnodes_std(femelement_table, ref_face_nodes)
|
|
for mf in ref_face_elements:
|
|
face_table[mf] = femelement_table[mf]
|
|
# print(face_table)
|
|
return face_table
|
|
|
|
|
|
def build_mesh_faces_of_volume_elements(face_table, femelement_table):
|
|
# node index of facenodes in femelementtable volume element
|
|
# if we know the position of the node we can build the element face out of the unsorted face nodes
|
|
face_nodenumber_table = {} # { volumeID : ( index, ... , index ) }
|
|
for veID in face_table:
|
|
face_nodenumber_table[veID] = []
|
|
for n in face_table[veID]:
|
|
index = femelement_table[veID].index(n)
|
|
# print(index)
|
|
face_nodenumber_table[veID].append(index + 1) # lokale node number = index + 1
|
|
# print('VolElement:', veID)
|
|
# print(' --> ', femelement_table[veID])
|
|
# print(' --> ', face_table[veID])
|
|
# print(' --> ', face_nodenumber_table[veID])
|
|
for veID in face_nodenumber_table:
|
|
vol_node_ct = len(femelement_table[veID])
|
|
face_node_indexs = sorted(face_nodenumber_table[veID])
|
|
if vol_node_ct == 10: # tetra10 --> tria6 face
|
|
if face_node_indexs == [1, 2, 3, 5, 6, 7]: # node order of face in tetra10 volume element
|
|
node_numbers = (1, 2, 3, 5, 6, 7) # node order of a tria6 face of tetra10
|
|
elif face_node_indexs == [1, 2, 4, 5, 8, 9]:
|
|
node_numbers = (1, 4, 2, 8, 9, 5)
|
|
elif face_node_indexs == [1, 3, 4, 7, 8, 10]:
|
|
node_numbers = (1, 3, 4, 7, 10, 8)
|
|
elif face_node_indexs == [2, 3, 4, 6, 9, 10]:
|
|
node_numbers = (2, 4, 3, 9, 10, 6)
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
|
|
elif vol_node_ct == 4: # tetra4 --> tria3 face
|
|
if face_node_indexs == [1, 2, 3]: # node order of face in tetra4 volume element
|
|
node_numbers = (1, 2, 3) # node order of a tria3 face of tetra4
|
|
elif face_node_indexs == [1, 2, 4]:
|
|
node_numbers = (1, 4, 2, 8)
|
|
elif face_node_indexs == [1, 3, 4]:
|
|
node_numbers = (1, 3, 4)
|
|
elif face_node_indexs == [2, 3, 4]:
|
|
node_numbers = (2, 4, 3)
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
|
|
elif vol_node_ct == 20: # hexa20 --> quad8 face
|
|
if face_node_indexs == [1, 2, 3, 4, 9, 10, 11, 12]: # node order of face in hexa20 volume element
|
|
node_numbers = (1, 2, 3, 4, 9, 10, 11, 12) # node order of a quad8 face of hexa20
|
|
elif face_node_indexs == [5, 6, 7, 8, 13, 14, 15, 16]:
|
|
node_numbers = (5, 8, 7, 6, 16, 15, 14, 13)
|
|
elif face_node_indexs == [1, 2, 5, 6, 9, 13, 17, 18]:
|
|
node_numbers = (1, 5, 6, 2, 17, 13, 18, 9)
|
|
elif face_node_indexs == [3, 4, 7, 8, 11, 15, 19, 20]:
|
|
node_numbers = (3, 7, 8, 4, 19, 15, 20, 11)
|
|
elif face_node_indexs == [1, 4, 5, 8, 12, 16, 17, 20]:
|
|
node_numbers = (1, 4, 8, 5, 12, 20, 16, 17)
|
|
elif face_node_indexs == [2, 3, 6, 7, 10, 14, 18, 19]:
|
|
node_numbers = (2, 6, 7, 3, 18, 14, 19, 10)
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
|
|
elif vol_node_ct == 8: # hexa8 --> quad4 face
|
|
face_node_indexs = sorted(face_nodenumber_table[veID])
|
|
if face_node_indexs == [1, 2, 3, 4]: # node order of face in hexa8 volume element
|
|
node_numbers = (1, 2, 3, 4) # node order of a quad8 face of hexa8
|
|
elif face_node_indexs == [5, 6, 7, 8]:
|
|
node_numbers = (5, 8, 7, 6)
|
|
elif face_node_indexs == [1, 2, 5, 6]:
|
|
node_numbers = (1, 5, 6, 2)
|
|
elif face_node_indexs == [3, 4, 7, 8]:
|
|
node_numbers = (3, 7, 8, 4)
|
|
elif face_node_indexs == [1, 4, 5, 8]:
|
|
node_numbers = (1, 4, 8, 5)
|
|
elif face_node_indexs == [2, 3, 6, 7]:
|
|
node_numbers = (2, 6, 7, 3)
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): hexa20: face not found!" + str(face_node_indexs) + "\n")
|
|
elif vol_node_ct == 15: # penta15 --> tria6 and quad8 faces
|
|
if face_node_indexs == [1, 2, 3, 7, 8, 9]: # node order of face in penta15 volume element
|
|
node_numbers = (1, 2, 3, 7, 8, 9) # node order of a tria6 face of penta15
|
|
elif face_node_indexs == [4, 5, 6, 10, 11, 12]:
|
|
node_numbers = (4, 6, 5, 12, 11, 10) # tria6
|
|
elif face_node_indexs == [1, 2, 4, 5, 7, 10, 13, 14]:
|
|
node_numbers = (1, 4, 5, 2, 13, 10, 14, 7) # quad8
|
|
elif face_node_indexs == [1, 3, 4, 6, 9, 12, 13, 15]:
|
|
node_numbers = (1, 3, 6, 4, 9, 15, 12, 13) # quad8
|
|
elif face_node_indexs == [2, 3, 5, 6, 8, 11, 14, 15]:
|
|
node_numbers = (2, 5, 6, 3, 14, 11, 15, 8) # quad8
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): penta15: face not found!" + str(face_node_indexs) + "\n")
|
|
elif vol_node_ct == 6: # penta6 --> tria3 and quad4 faces
|
|
if face_node_indexs == [1, 2, 3]: # node order of face in penta6 volume element
|
|
node_numbers = (1, 2, 3) # node order of a tria3 face of penta6
|
|
elif face_node_indexs == [4, 5, 6]:
|
|
node_numbers = (4, 6, 5) # tria3
|
|
elif face_node_indexs == [1, 2, 4, 5]:
|
|
node_numbers = (1, 4, 5, 2) # quad4
|
|
elif face_node_indexs == [1, 3, 4, 6]:
|
|
node_numbers = (1, 3, 6, 4) # quad4
|
|
elif face_node_indexs == [2, 3, 5, 6]:
|
|
node_numbers = (2, 5, 6, 3) # quad4
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): pent6: face not found!" + str(face_node_indexs) + "\n")
|
|
else:
|
|
FreeCAD.Console.PrintError("Error in build_mesh_faces_of_volume_elements(): Volume not implemented: volume node count" + str(vol_node_ct) + "\n")
|
|
face_nodes = []
|
|
for i in node_numbers:
|
|
i -= 1 # node_number starts with 1, index starts with 0 --> index = node number - 1
|
|
face_nodes.append(femelement_table[veID][i])
|
|
face_table[veID] = face_nodes # reset the entry in face_table
|
|
# print(' --> ', face_table[veID])
|
|
return face_table
|
|
|
|
|
|
def get_ref_facenodes_areas(femnodes_mesh, face_table):
|
|
# calulate the appropriate node_areas for every node of every mesh face (mf)
|
|
# G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
|
|
# FIXME only gives exact results in case of a real triangle. If for S6 or C3D10 elements
|
|
# the midnodes are not on the line between the end nodes the area will not be a triangle
|
|
# see http://forum.freecadweb.org/viewtopic.php?f=18&t=10939&start=40#p91355 and ff
|
|
# same applies for the quads, results are exact only if mid nodes are on the line between corner nodes
|
|
|
|
# [ (nodeID,Area), ... , (nodeID,Area) ] some nodes will have more than one entry
|
|
if (not femnodes_mesh) or (not face_table):
|
|
FreeCAD.Console.PrintError("Error: Empty femnodes_mesh or face_table!\n")
|
|
return []
|
|
node_area_table = []
|
|
mesh_face_area = 0
|
|
for mf in face_table:
|
|
femmesh_facetype = len(face_table[mf])
|
|
# nodes in face_table need to be in the right node order for the following calcualtions
|
|
if femmesh_facetype == 3: # 3 node femmesh face triangle
|
|
# corner_node_area = mesh_face_area / 3.0
|
|
# P3
|
|
# /\
|
|
# / \
|
|
# /____\
|
|
# P1 P2
|
|
P1 = femnodes_mesh[face_table[mf][0]]
|
|
P2 = femnodes_mesh[face_table[mf][1]]
|
|
P3 = femnodes_mesh[face_table[mf][2]]
|
|
|
|
mesh_face_area = get_triangle_area(P1, P2, P3)
|
|
corner_node_area = mesh_face_area / 3.0
|
|
|
|
node_area_table.append((face_table[mf][0], corner_node_area))
|
|
node_area_table.append((face_table[mf][1], corner_node_area))
|
|
node_area_table.append((face_table[mf][2], corner_node_area))
|
|
|
|
elif femmesh_facetype == 4: # 4 node femmesh face quad
|
|
# corner_node_area = mesh_face_area / 4.0
|
|
# P4_______P3
|
|
# | /|
|
|
# | t2 / |
|
|
# | / |
|
|
# | / |
|
|
# | / t1 |
|
|
# |/_____|
|
|
# P1 P2
|
|
P1 = femnodes_mesh[face_table[mf][0]]
|
|
P2 = femnodes_mesh[face_table[mf][1]]
|
|
P3 = femnodes_mesh[face_table[mf][2]]
|
|
P4 = femnodes_mesh[face_table[mf][3]]
|
|
|
|
mesh_face_t1_area = get_triangle_area(P1, P2, P3)
|
|
mesh_face_t2_area = get_triangle_area(P1, P3, P4)
|
|
mesh_face_area = mesh_face_t1_area + mesh_face_t1_area
|
|
corner_node_area = mesh_face_area / 4.0
|
|
|
|
node_area_table.append((face_table[mf][0], corner_node_area))
|
|
node_area_table.append((face_table[mf][1], corner_node_area))
|
|
node_area_table.append((face_table[mf][2], corner_node_area))
|
|
node_area_table.append((face_table[mf][3], corner_node_area))
|
|
|
|
elif femmesh_facetype == 6: # 6 node femmesh face triangle
|
|
# corner_node_area = 0
|
|
# middle_node_area = mesh_face_area / 3.0
|
|
# P3
|
|
# /\
|
|
# /t3\
|
|
# / \
|
|
# P6------P5
|
|
# / \ t4 / \
|
|
# /t1 \ /t2 \
|
|
# /_____\/_____\
|
|
# P1 P4 P2
|
|
P1 = femnodes_mesh[face_table[mf][0]]
|
|
P2 = femnodes_mesh[face_table[mf][1]]
|
|
P3 = femnodes_mesh[face_table[mf][2]]
|
|
P4 = femnodes_mesh[face_table[mf][3]]
|
|
P5 = femnodes_mesh[face_table[mf][4]]
|
|
P6 = femnodes_mesh[face_table[mf][5]]
|
|
|
|
mesh_face_t1_area = get_triangle_area(P1, P4, P6)
|
|
mesh_face_t2_area = get_triangle_area(P2, P5, P4)
|
|
mesh_face_t3_area = get_triangle_area(P3, P6, P5)
|
|
mesh_face_t4_area = get_triangle_area(P4, P5, P6)
|
|
mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area
|
|
middle_node_area = mesh_face_area / 3.0
|
|
|
|
node_area_table.append((face_table[mf][0], 0))
|
|
node_area_table.append((face_table[mf][1], 0))
|
|
node_area_table.append((face_table[mf][2], 0))
|
|
node_area_table.append((face_table[mf][3], middle_node_area))
|
|
node_area_table.append((face_table[mf][4], middle_node_area))
|
|
node_area_table.append((face_table[mf][5], middle_node_area))
|
|
|
|
elif femmesh_facetype == 8: # 8 node femmesh face quad
|
|
# corner_node_area = -mesh_face_area / 12.0 (negativ!)
|
|
# mid-side nodes = mesh_face_area / 3.0
|
|
# P4_________P7________P3
|
|
# | / | \ |
|
|
# | t4 / | \ t3 |
|
|
# | / | \ |
|
|
# |/ | \|
|
|
# P8| t5 | t6 |P6
|
|
# |\ | / |
|
|
# | \ | / |
|
|
# | t1 \ | / t2 |
|
|
# |______\_|_/_______|
|
|
# P1 P5 P2
|
|
P1 = femnodes_mesh[face_table[mf][0]]
|
|
P2 = femnodes_mesh[face_table[mf][1]]
|
|
P3 = femnodes_mesh[face_table[mf][2]]
|
|
P4 = femnodes_mesh[face_table[mf][3]]
|
|
P5 = femnodes_mesh[face_table[mf][4]]
|
|
P6 = femnodes_mesh[face_table[mf][5]]
|
|
P7 = femnodes_mesh[face_table[mf][6]]
|
|
P8 = femnodes_mesh[face_table[mf][7]]
|
|
|
|
mesh_face_t1_area = get_triangle_area(P1, P5, P8)
|
|
mesh_face_t2_area = get_triangle_area(P5, P2, P6)
|
|
mesh_face_t3_area = get_triangle_area(P6, P3, P7)
|
|
mesh_face_t4_area = get_triangle_area(P7, P4, P8)
|
|
mesh_face_t5_area = get_triangle_area(P5, P7, P8)
|
|
mesh_face_t6_area = get_triangle_area(P5, P6, P7)
|
|
mesh_face_area = mesh_face_t1_area + mesh_face_t2_area + mesh_face_t3_area + mesh_face_t4_area + mesh_face_t5_area + mesh_face_t6_area
|
|
corner_node_area = -mesh_face_area / 12.0
|
|
middle_node_area = mesh_face_area / 3.0
|
|
|
|
node_area_table.append((face_table[mf][0], corner_node_area))
|
|
node_area_table.append((face_table[mf][1], corner_node_area))
|
|
node_area_table.append((face_table[mf][2], corner_node_area))
|
|
node_area_table.append((face_table[mf][3], corner_node_area))
|
|
node_area_table.append((face_table[mf][4], middle_node_area))
|
|
node_area_table.append((face_table[mf][5], middle_node_area))
|
|
node_area_table.append((face_table[mf][6], middle_node_area))
|
|
node_area_table.append((face_table[mf][7], middle_node_area))
|
|
return node_area_table
|
|
|
|
|
|
def get_ref_shape_node_sum_geom_table(node_geom_table):
|
|
# shape could be Edge or Face, geom could be lenght or area
|
|
# summ of legth or area for each node of the ref_shape
|
|
node_sum_geom_table = {}
|
|
for n, A in node_geom_table:
|
|
# print(n, ' --> ', A)
|
|
if n in node_sum_geom_table:
|
|
node_sum_geom_table[n] = node_sum_geom_table[n] + A
|
|
else:
|
|
node_sum_geom_table[n] = A
|
|
return node_sum_geom_table
|
|
|
|
|
|
def get_analysis_group_elements(aAnalysis, aPart):
|
|
''' all Reference shapes of all Analysis member are searched in the Shape of aPart. If found in shape they are added to a dict
|
|
{ConstraintName : ['ShapeType of the Elements'], [ElementID, ElementID, ...], ...}
|
|
'''
|
|
aShape = aPart.Shape
|
|
group_elements = {} # { name : [element, element, ... , element]}
|
|
empty_references = []
|
|
for m in aAnalysis.Member:
|
|
if hasattr(m, "References"):
|
|
# print(m.Name)
|
|
key = m.Name
|
|
elements = []
|
|
stype = None
|
|
if m.References:
|
|
for r in m.References:
|
|
parent = r[0]
|
|
childs = r[1]
|
|
# print(parent)
|
|
# print(childs)
|
|
for child in childs:
|
|
ref_shape = get_element(parent, child) # the method getElement(element) does not return Solid elements
|
|
if not stype:
|
|
stype = ref_shape.ShapeType
|
|
elif stype != ref_shape.ShapeType:
|
|
FreeCAD.Console.PrintError('Error, two refschapes in References with different ShapeTypes.\n')
|
|
# print(ref_shape)
|
|
found_element = find_element_in_shape(aShape, ref_shape)
|
|
if found_element is not None:
|
|
elements.append(found_element)
|
|
else:
|
|
FreeCAD.Console.PrintError('Problem: No element found for: ' + str(ref_shape) + '\n')
|
|
print(' ' + m.Name)
|
|
print(' ' + str(m.References))
|
|
print(' ' + r[0].Name)
|
|
group_elements[key] = sorted(elements)
|
|
else:
|
|
print(' Empty reference: ' + m.Name)
|
|
empty_references.append(m)
|
|
if empty_references:
|
|
if len(empty_references) == 1:
|
|
group_elements = get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShape)
|
|
else:
|
|
FreeCAD.Console.PrintError('Problem: more than one object with empty references.\n')
|
|
print('We gone try to get the empty material references anyway.\n')
|
|
# ShellThickness and BeamSection could have empty references, but on solid meshes only materials should have empty references
|
|
for er in empty_references:
|
|
print(er.Name)
|
|
group_elements = get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShape)
|
|
# check if all groups have elements:
|
|
for g in group_elements:
|
|
# print(group_elements[g])
|
|
if len(group_elements[g]) == 0:
|
|
FreeCAD.Console.PrintError('Error: shapes for: ' + g + 'not found!\n')
|
|
return group_elements
|
|
|
|
|
|
def get_anlysis_empty_references_group_elements(group_elements, aAnalysis, aShape):
|
|
'''get the elementIDs if the Reference shape is empty
|
|
see get_analysis_group_elements() for more informatations
|
|
on solid meshes only material objects could have an empty reference without beeing something wrong!
|
|
face meshes could have empty ShellThickness and edge meshes could have empty BeamSection
|
|
'''
|
|
# print(group_elements)
|
|
material_ref_shapes = []
|
|
material_shape_type = ''
|
|
missed_material_refshapes = []
|
|
empty_reference_material = None
|
|
for m in aAnalysis.Member:
|
|
if m.isDerivedFrom("App::MaterialObjectPython"):
|
|
if hasattr(m, "References") and not m.References:
|
|
if not empty_reference_material:
|
|
empty_reference_material = m.Name
|
|
else:
|
|
FreeCAD.Console.PrintError('Problem in get_anlysis_empty_references_group_elements, we seams to have two or more materials with empty referneces')
|
|
return {}
|
|
elif hasattr(m, "References") and m.References:
|
|
# ShapeType ot the group elements, strip the number of the first group element
|
|
# http://stackoverflow.com/questions/12851791/removing-numbers-from-string
|
|
group_shape_type = ''.join(i for i in group_elements[m.Name][0] if not i.isdigit())
|
|
if not material_shape_type:
|
|
material_shape_type = group_shape_type
|
|
elif material_shape_type != group_shape_type:
|
|
FreeCAD.Console.PrintError('Problem, material shape type does not match get_anlysis_empty_references_group_elements')
|
|
for ele in group_elements[m.Name]:
|
|
material_ref_shapes.append(ele)
|
|
if material_shape_type == 'Solid':
|
|
# print(len(aShape.Solids))
|
|
for i in range(len(aShape.Solids)):
|
|
ele = 'Solid' + str(i + 1)
|
|
if ele not in material_ref_shapes:
|
|
missed_material_refshapes.append(ele)
|
|
elif material_shape_type == 'Face':
|
|
# print(len(aShape.Faces))
|
|
for i in range(len(aShape.Faces)):
|
|
ele = 'Face' + str(i + 1)
|
|
if ele not in material_ref_shapes:
|
|
missed_material_refshapes.append(ele)
|
|
elif material_shape_type == 'Edge':
|
|
# print(len(aShape.Edges))
|
|
for i in range(len(aShape.Edges)):
|
|
ele = 'Edge' + str(i + 1)
|
|
if ele not in material_ref_shapes:
|
|
missed_material_refshapes.append(ele)
|
|
else:
|
|
print(' One material with no reference shapes. No need to make a group for materials.')
|
|
# make no changes group_elements
|
|
return group_elements
|
|
# print(sorted(material_ref_shapes))
|
|
# print(sorted(missed_material_refshapes))
|
|
# print(group_elements)
|
|
group_elements[empty_reference_material] = sorted(missed_material_refshapes)
|
|
# print(group_elements)
|
|
return group_elements
|
|
|
|
|
|
def find_element_in_shape(aShape, anElement):
|
|
# import Part
|
|
ele_st = anElement.ShapeType
|
|
if ele_st == 'Solid' or ele_st == 'CompSolid':
|
|
for index, solid in enumerate(aShape.Solids):
|
|
# print(is_same_geometry(solid, anElement))
|
|
if is_same_geometry(solid, anElement):
|
|
# print(index)
|
|
# Part.show(aShape.Solids[index])
|
|
ele = ele_st + str(index + 1)
|
|
return ele
|
|
FreeCAD.Console.PrintError('Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n')
|
|
if ele_st == 'Solid' and aShape.ShapeType == 'Solid':
|
|
print('We have been searching for a Solid in a Solid and we have not found it. In most cases this should be searching for a Solid inside a CompSolid. Check the ShapeType of your Part to mesh.')
|
|
# Part.show(anElement)
|
|
# Part.show(aShape)
|
|
elif ele_st == 'Face' or ele_st == 'Shell':
|
|
for index, face in enumerate(aShape.Faces):
|
|
# print(is_same_geometry(face, anElement))
|
|
if is_same_geometry(face, anElement):
|
|
# print(index)
|
|
# Part.show(aShape.Faces[index])
|
|
ele = ele_st + str(index + 1)
|
|
return ele
|
|
elif ele_st == 'Edge' or ele_st == 'Wire':
|
|
for index, edge in enumerate(aShape.Edges):
|
|
# print(is_same_geometry(edge, anElement))
|
|
if is_same_geometry(edge, anElement):
|
|
# print(index)
|
|
# Part.show(aShape.Edges[index])
|
|
ele = ele_st + str(index + 1)
|
|
return ele
|
|
elif ele_st == 'Vertex':
|
|
for index, vertex in enumerate(aShape.Vertexes):
|
|
# print(is_same_geometry(vertex, anElement))
|
|
if is_same_geometry(vertex, anElement):
|
|
# print(index)
|
|
# Part.show(aShape.Vertexes[index])
|
|
ele = ele_st + str(index + 1)
|
|
return ele
|
|
elif ele_st == 'Compound':
|
|
FreeCAD.Console.PrintError('Compound is not supported.\n')
|
|
|
|
|
|
def get_vertexes_by_element(aShape, anElement):
|
|
# we gone extent the method find_element_in_shape and return the vertexes
|
|
# import Part
|
|
ele_vertexes = []
|
|
ele_st = anElement.ShapeType
|
|
if ele_st == 'Solid' or ele_st == 'CompSolid':
|
|
for index, solid in enumerate(aShape.Solids):
|
|
if is_same_geometry(solid, anElement):
|
|
for vele in aShape.Solids[index].Vertexes:
|
|
for i, v in enumerate(aShape.Vertexes):
|
|
if vele.isSame(v): # use isSame, because orientation could be different
|
|
ele_vertexes.append(i)
|
|
# print(' ' + str(sorted(ele_vertexes)))
|
|
return ele_vertexes
|
|
FreeCAD.Console.PrintError('Error, Solid ' + str(anElement) + ' not found in: ' + str(aShape) + '\n')
|
|
elif ele_st == 'Face' or ele_st == 'Shell':
|
|
for index, face in enumerate(aShape.Faces):
|
|
if is_same_geometry(face, anElement):
|
|
for vele in aShape.Faces[index].Vertexes:
|
|
for i, v in enumerate(aShape.Vertexes):
|
|
if vele.isSame(v): # use isSame, because orientation could be different
|
|
ele_vertexes.append(i)
|
|
# print(' ' + str(sorted(ele_vertexes)))
|
|
return ele_vertexes
|
|
elif ele_st == 'Edge' or ele_st == 'Wire':
|
|
for index, edge in enumerate(aShape.Edges):
|
|
if is_same_geometry(edge, anElement):
|
|
for vele in aShape.Edges[index].Vertexes:
|
|
for i, v in enumerate(aShape.Vertexes):
|
|
if vele.isSame(v): # use isSame, because orientation could be different
|
|
ele_vertexes.append(i)
|
|
# print(' ' + str(sorted(ele_vertexes)))
|
|
return ele_vertexes
|
|
elif ele_st == 'Vertex':
|
|
for index, vertex in enumerate(aShape.Vertexes):
|
|
if is_same_geometry(vertex, anElement):
|
|
ele_vertexes.append(index)
|
|
# print(' ' + str(sorted(ele_vertexes)))
|
|
return ele_vertexes
|
|
elif ele_st == 'Compound':
|
|
FreeCAD.Console.PrintError('Compound is not supported.\n')
|
|
|
|
|
|
def is_same_geometry(shape1, shape2):
|
|
# the vertexes and the CenterOfMass are compared
|
|
# it is a hack, but I do not know any better !
|
|
# check of Volume and Area before starting with the vertices could be added
|
|
# BoundBox is possible too, but is BB calcualtions robust?!
|
|
# print(shape1)
|
|
# print(shape2)
|
|
same_Vertexes = 0
|
|
if len(shape1.Vertexes) == len(shape2.Vertexes) and len(shape1.Vertexes) > 1:
|
|
# compare CenterOfMass
|
|
if shape1.CenterOfMass != shape2.CenterOfMass:
|
|
return False
|
|
else:
|
|
# compare the Vertexes
|
|
for vs1 in shape1.Vertexes:
|
|
for vs2 in shape2.Vertexes:
|
|
if vs1.X == vs2.X and vs1.Y == vs2.Y and vs1.Z == vs2.Z:
|
|
same_Vertexes += 1
|
|
continue
|
|
# print(same_Vertexes)
|
|
if same_Vertexes == len(shape1.Vertexes):
|
|
return True
|
|
else:
|
|
return False
|
|
if len(shape1.Vertexes) == len(shape2.Vertexes) and len(shape1.Vertexes) == 1:
|
|
vs1 = shape1.Vertexes[0]
|
|
vs2 = shape2.Vertexes[0]
|
|
if vs1.X == vs2.X and vs1.Y == vs2.Y and vs1.Z == vs2.Z:
|
|
return True
|
|
else:
|
|
return False
|
|
else:
|
|
return False
|
|
|
|
|
|
def get_element(part, element):
|
|
if element.startswith('Solid'):
|
|
return part.Shape.Solids[int(element.lstrip('Solid')) - 1] # Solid
|
|
else:
|
|
return part.Shape.getElement(element) # Face, Edge, Vertex
|
|
|
|
|
|
def femelements_count_ok(len_femelement_table, count_femelements):
|
|
if count_femelements == len_femelement_table:
|
|
print('Count FEM elements as sum of constraints: ', count_femelements)
|
|
print('Count FEM elements of the FreeCAD FEM mesh: ', len_femelement_table)
|
|
return True
|
|
else:
|
|
print('ERROR: femelement_table != count_femelements')
|
|
print('Count FEM elements as sum of constraints: ', count_femelements)
|
|
print('Count FEM elements of the FreeCAD FEM Mesh: ', len_femelement_table)
|
|
return False
|
|
|
|
|
|
def delete_duplicate_mesh_elements(refelement_table):
|
|
new_refelement_table = {} # duplicates deleted
|
|
for elem, nodes in refelement_table.items():
|
|
if sorted(nodes) not in sortlistoflistvalues(new_refelement_table.values()):
|
|
new_refelement_table[elem] = nodes
|
|
return new_refelement_table
|
|
|
|
|
|
def get_triangle_area(P1, P2, P3):
|
|
# import Part
|
|
# W = Part.Wire([Part.makeLine(P1,P2), Part.makeLine(P2,P3), Part.makeLine(P3,P1)])
|
|
# Part.show(Part.Face(W))
|
|
vec1 = P2 - P1
|
|
vec2 = P3 - P1
|
|
vec3 = vec1.cross(vec2)
|
|
return 0.5 * vec3.Length
|
|
|
|
|
|
def sortlistoflistvalues(listoflists):
|
|
new_list = []
|
|
for l in listoflists:
|
|
new_list.append(sorted(l))
|
|
return new_list
|
|
|
|
|
|
def is_solid_femmesh(femmesh):
|
|
if femmesh.VolumeCount > 0: # solid femmesh
|
|
return True
|
|
|
|
|
|
def has_no_face_data(femmesh):
|
|
if femmesh.FaceCount == 0: # femmesh has no face data, could be a edge femmesh or a solid femmesh without face data
|
|
return True
|
|
|
|
|
|
def is_face_femmesh(femmesh):
|
|
if femmesh.VolumeCount == 0 and femmesh.FaceCount > 0: # face femmesh
|
|
return True
|
|
|
|
|
|
def is_edge_femmesh(femmesh):
|
|
if femmesh.VolumeCount == 0 and femmesh.FaceCount == 0 and femmesh.EdgeCount > 0: # edge femmesh
|
|
return True
|
|
|
|
|
|
def get_three_non_colinear_nodes(nodes_coords):
|
|
# Code to obtain three non-colinear nodes on the PlaneRotation support face
|
|
# nodes_coords --> [(nodenumber, x, y, z), (nodenumber, x, y, z), ...]
|
|
if not nodes_coords:
|
|
print(len(nodes_coords))
|
|
print('Error: No nodes in nodes_coords')
|
|
return []
|
|
dum_max = [1, 2, 3, 4, 5, 6, 7, 8, 0]
|
|
for i in range(len(nodes_coords)):
|
|
for j in range(len(nodes_coords) - 1 - i):
|
|
x_1 = nodes_coords[j][1]
|
|
x_2 = nodes_coords[j + 1][1]
|
|
y_1 = nodes_coords[j][2]
|
|
y_2 = nodes_coords[j + 1][2]
|
|
z_1 = nodes_coords[j][3]
|
|
z_2 = nodes_coords[j + 1][3]
|
|
node_1 = nodes_coords[j][0]
|
|
node_2 = nodes_coords[j + 1][0]
|
|
distance = ((x_1 - x_2) ** 2 + (y_1 - y_2) ** 2 + (z_1 - z_2) ** 2) ** 0.5
|
|
if distance > dum_max[8]:
|
|
dum_max = [node_1, x_1, y_1, z_1, node_2, x_2, y_2, z_2, distance]
|
|
node_dis = [1, 0]
|
|
for i in range(len(nodes_coords)):
|
|
x_1 = dum_max[1]
|
|
x_2 = dum_max[5]
|
|
x_3 = nodes_coords[i][1]
|
|
y_1 = dum_max[2]
|
|
y_2 = dum_max[6]
|
|
y_3 = nodes_coords[i][2]
|
|
z_1 = dum_max[3]
|
|
z_2 = dum_max[7]
|
|
z_3 = nodes_coords[i][3]
|
|
node_3 = int(nodes_coords[j][0])
|
|
distance_1 = ((x_1 - x_3) ** 2 + (y_1 - y_3) ** 2 + (z_1 - z_3) ** 2) ** 0.5
|
|
distance_2 = ((x_3 - x_2) ** 2 + (y_3 - y_2) ** 2 + (z_3 - z_2) ** 2) ** 0.5
|
|
tot = distance_1 + distance_2
|
|
if tot > node_dis[1]:
|
|
node_dis = [node_3, tot]
|
|
node_1 = int(dum_max[0])
|
|
node_2 = int(dum_max[4])
|
|
print([node_1, node_2, node_3])
|
|
return [node_1, node_2, node_3]
|
|
|
|
|
|
def get_rectangular_coords(obj):
|
|
from math import cos, sin, radians
|
|
A = [1, 0, 0]
|
|
B = [0, 1, 0]
|
|
a_x = A[0]
|
|
a_y = A[1]
|
|
a_z = A[2]
|
|
b_x = A[0]
|
|
b_y = A[1]
|
|
b_z = A[2]
|
|
x_rot = radians(obj.X_rot)
|
|
y_rot = radians(obj.Y_rot)
|
|
z_rot = radians(obj.Z_rot)
|
|
if obj.X_rot != 0:
|
|
a_x = A[0]
|
|
a_y = A[1] * cos(x_rot) + A[2] * sin(x_rot)
|
|
a_z = A[2] * cos(x_rot) - A[1] * sin(x_rot)
|
|
b_x = B[0]
|
|
b_y = B[1] * cos(x_rot) + B[2] * sin(x_rot)
|
|
b_z = B[2] * cos(x_rot) - B[1] * sin(x_rot)
|
|
A = [a_x, a_y, a_z]
|
|
B = [b_x, b_y, b_z]
|
|
if obj.Y_rot != 0:
|
|
a_x = A[0] * cos(y_rot) - A[2] * sin(y_rot)
|
|
a_y = A[1]
|
|
a_z = A[2] * cos(y_rot) + A[0] * sin(y_rot)
|
|
b_x = B[0] * cos(y_rot) - B[2] * sin(y_rot)
|
|
b_y = B[1]
|
|
b_z = B[2] * cos(y_rot) + B[0] * sin(z_rot)
|
|
A = [a_x, a_y, a_z]
|
|
B = [b_x, b_y, b_z]
|
|
if obj.Z_rot != 0:
|
|
a_x = A[0] * cos(z_rot) + A[1] * sin(z_rot)
|
|
a_y = A[1] * cos(z_rot) - A[0] * sin(z_rot)
|
|
a_z = A[2]
|
|
b_x = B[0] * cos(z_rot) + B[1] * sin(z_rot)
|
|
b_y = B[1] * cos(z_rot) - B[0] * sin(z_rot)
|
|
b_z = B[2]
|
|
A = [a_x, a_y, a_z]
|
|
B = [b_x, b_y, b_z]
|
|
|
|
A_coords = str(round(A[0], 4)) + ',' + str(round(A[1], 4)) + ',' + str(round(A[2], 4))
|
|
B_coords = str(round(B[0], 4)) + ',' + str(round(B[1], 4)) + ',' + str(round(B[2], 4))
|
|
coords = A_coords + ',' + B_coords
|
|
return coords
|
|
|
|
|
|
def get_cylindrical_coords(obj):
|
|
vec = obj.Axis
|
|
base = obj.BasePoint
|
|
Ax = base[0] + 10 * vec[0]
|
|
Ay = base[1] + 10 * vec[1]
|
|
Az = base[2] + 10 * vec[2]
|
|
Bx = base[0] - 10 * vec[0]
|
|
By = base[1] - 10 * vec[1]
|
|
Bz = base[2] - 10 * vec[2]
|
|
A = [Ax, Ay, Az]
|
|
B = [Bx, By, Bz]
|
|
A_coords = str(A[0]) + ',' + str(A[1]) + ',' + str(A[2])
|
|
B_coords = str(B[0]) + ',' + str(B[1]) + ',' + str(B[2])
|
|
coords = A_coords + ',' + B_coords
|
|
return coords
|
|
|
|
|
|
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
|
|
('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])
|
|
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".format(
|
|
len(elms_hexa20), len(elms_tria3), len(elms_tria6), len(elms_quad4), len(elms_quad8), len(elms_seg2)))
|
|
else:
|
|
FreeCAD.Console.PrintError("No Elements found!\n")
|
|
else:
|
|
FreeCAD.Console.PrintError("No Nodes found!\n")
|
|
return mesh
|
|
|
|
# @}
|