Skip to content

Commit

Permalink
FEM: ccxwriter, much more exact results for cload on edges
Browse files Browse the repository at this point in the history
  • Loading branch information
@berndhahnebach @wwmayer
berndhahnebach authored and wwmayer committed Oct 26, 2015
1 parent 9c0bcb7 commit e5adec9
Showing 1 changed file with 157 additions and 21 deletions.
178 changes: 157 additions & 21 deletions src/Mod/Fem/ccxInpWriter.py
View file
Expand Up @@ -167,27 +167,24 @@ def write_node_sets_constraints_force(self, f):
f.write('** Node sets for loads\n')
f.write('** written by {} function\n'.format(sys._getframe().f_code.co_name))
for fobj in self.force_objects:
frc_obj = fobj['Object']
f.write('*NSET,NSET=' + frc_obj.Name + '\n')
NbrForceNodes = 0
for o, elem in frc_obj.References:
fo = o.Shape.getElement(elem)
n = []
if fobj['RefShapeType'] == 'Vertex':
if fobj['RefShapeType'] == 'Vertex':
frc_obj = fobj['Object']
f.write('*NSET,NSET=' + frc_obj.Name + '\n')
NbrForceNodes = 0
for o, elem in frc_obj.References:
fo = o.Shape.getElement(elem)
n = []
n = self.mesh_object.FemMesh.getNodesByVertex(fo)
elif fobj['RefShapeType'] == 'Edge':
n = self.mesh_object.FemMesh.getNodesByEdge(fo)
for i in n:
f.write(str(i) + ',\n')
NbrForceNodes = NbrForceNodes + 1 # NodeSum of mesh-nodes of ALL reference shapes from force_object
# calculate node load
if NbrForceNodes != 0:
fobj['NodeLoad'] = (frc_obj.Force) / NbrForceNodes
# FIXME for loads on edges the node count is used to distribute the load on the edges.
# In case of a not uniform fem mesh this could result in wrong force distribution
# and thus in wrong analysis results. see def write_constraints_force()
f.write('** concentrated load [N] distributed on all mesh nodes of the given shapes\n')
f.write('** ' + str(frc_obj.Force) + ' N / ' + str(NbrForceNodes) + ' Nodes = ' + str(fobj['NodeLoad']) + ' N on each node\n')
for i in n:
f.write(str(i) + ',\n')
NbrForceNodes = NbrForceNodes + 1 # NodeSum of mesh-nodes of ALL reference shapes from force_object
# calculate node load
if NbrForceNodes != 0:
fobj['NodeLoad'] = (frc_obj.Force) / NbrForceNodes
f.write('** concentrated load [N] distributed on all mesh nodes of the given vertieces\n')
f.write('** ' + str(frc_obj.Force) + ' N / ' + str(NbrForceNodes) + ' Nodes = ' + str(fobj['NodeLoad']) + ' N on each node\n')
else:
f.write('** no point load on vertices --> no set for node loads\n')

def write_materials(self, f):
f.write('\n***********************************************************\n')
Expand Down Expand Up @@ -274,7 +271,7 @@ def write_constraints_force(self, f):
if frc_obj.Force == 0:
print(' Warning --> Force = 0')

if fobj['RefShapeType'] == 'Vertex' or fobj['RefShapeType'] == 'Edge': # load on edges or vertieces
if fobj['RefShapeType'] == 'Vertex': # point load on vertieces
node_load = fobj['NodeLoad']
frc_obj_name = frc_obj.Name
f.write('** force: ' + str(node_load) + ' N, direction: ' + str(direction_vec) + '\n')
Expand All @@ -285,6 +282,61 @@ def write_constraints_force(self, f):
f.write(frc_obj_name + ',2,' + v2 + '\n')
f.write(frc_obj_name + ',3,' + v3 + '\n\n')

elif fobj['RefShapeType'] == 'Edge': # line load on edges
sum_ref_edge_length = 0
sum_ref_edge_node_length = 0 # for debugging
sum_node_load = 0 # for debugging
for o, elem in frc_obj.References:
elem_o = o.Shape.getElement(elem)
sum_ref_edge_length += elem_o.Length
if sum_ref_edge_length != 0:
force_per_sum_ref_edge_length = frc_obj.Force / sum_ref_edge_length
for o, elem in frc_obj.References:
elem_o = o.Shape.getElement(elem)
ref_edge = elem_o

# edge_table = { meshedgeID : ( nodeID, ... , nodeID ) }
edge_table = self.get_refedge_node_table(ref_edge)

# node_length_table = [ (nodeID, length), ... , (nodeID, length) ] some nodes will have more than one entry
node_length_table = self.get_refedge_node_lengths(edge_table)

# node_sum_length_table = { nodeID : Length, ... , nodeID : Length } LengthSum for each node, one entry for each node
node_sum_length_table = self.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
sum_ref_edge_node_length += sum_node_lengths

f.write('** node loads on element ' + fobj['RefShapeType'] + ': ' + o.Name + ':' + elem + '\n')
for n in sorted(node_load_table):
node_load = node_load_table[n]
sum_node_load += node_load # for debugging
if (direction_vec.x != 0.0):
v1 = "{:.13E}".format(direction_vec.x * node_load)
f.write(str(n) + ',1,' + v1 + '\n')
if (direction_vec.y != 0.0):
v2 = "{:.13E}".format(direction_vec.y * node_load)
f.write(str(n) + ',2,' + v2 + '\n')
if (direction_vec.z != 0.0):
v3 = "{:.13E}".format(direction_vec.z * node_load)
f.write(str(n) + ',3,' + v3 + '\n')
f.write('\n')
f.write('\n')
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')

elif fobj['RefShapeType'] == 'Face': # area load on faces
sum_ref_face_area = 0
sum_ref_face_node_area = 0 # for debugging
Expand Down Expand Up @@ -594,6 +646,49 @@ def get_femelement_sets(self, fem_objects):
if not femelements_count_ok(self.fem_element_table, count_femelements):
FreeCAD.Console.PrintError('Error in get_femelement_sets -- > femelements_count_ok failed!\n')

def get_refedge_node_table(self, refedge):
edge_table = {} # { meshedgeID : ( nodeID, ... , nodeID ) }
refedge_nodes = self.mesh_object.FemMesh.getNodesByEdge(refedge)
if is_solid_mesh(self.mesh_object.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 self.fem_element_table:
nodecount = 0
for node in self.fem_element_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 self.fem_element_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
elif is_shell_mesh(self.mesh_object.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 self.fem_element_table:
nodecount = 0
for node in self.fem_element_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 self.fem_element_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
elif is_beam_mesh(self.mesh_object.FemMesh):
refedge_fem_edgeelements = getFemElementsByNodes(self.fem_element_table, refedge_nodes)
for elem in refedge_fem_edgeelements:
edge_table[elem] = self.fem_element_table[elem] # { edgeID : ( nodeID, ... , nodeID )} # all nodes off this femedgeelement
return edge_table

def get_ref_face_node_table(self, ref_face):
face_table = {} # { meshfaceID : ( nodeID, ... , nodeID ) }
if is_solid_mesh(self.mesh_object.FemMesh):
Expand All @@ -619,6 +714,47 @@ def get_ref_face_node_table(self, ref_face):
face_table[mf] = self.fem_element_table[mf]
return face_table

def get_refedge_node_lengths(self, 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
node_length_table = []
mesh_edge_length = 0
# print len(edge_table)
for me in edge_table:
if len(edge_table[me]) == 2: # 2 node mesh edge
# end_node_length = mesh_edge_length / 2
# ______
# P1 P2
P1 = self.mesh_object.FemMesh.Nodes[edge_table[me][0]]
P2 = self.mesh_object.FemMesh.Nodes[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 len(edge_table[me]) == 3: # 3 node mesh edge
# end_node_length = mesh_edge_length / 6
# middle_node_length = mesh_face_area * 2 / 3
# _______ _______
# P1 P3 P2
P1 = self.mesh_object.FemMesh.Nodes[edge_table[me][0]]
P2 = self.mesh_object.FemMesh.Nodes[edge_table[me][1]]
P3 = self.mesh_object.FemMesh.Nodes[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_face_node_areas(self, face_table):
# calulate the appropriate node_areas for every node of every mesh face (mf)
# G. Lakshmi Narasaiah, Finite Element Analysis, p206ff
Expand Down

0 comments on commit e5adec9

Please sign in to comment.