mr_sceptical athotmail.com
2011-05-12 15:35:16 UTC
Dear All
I am trying to write a Python script to assign beam section orientation to a random network of beams. The part is made up of individual beams (in 3d space) and merged into a single part. Intersections are retained so that different parts can be assigned different properties. This causes an issue, as to assign beam section orientation to the part (necessary once it has been merged), it is necessary to select all of the regions between the intersections. If a single beam is intersected by one or more other beams, each part must be selected and an orientation assigned.
My current process for this is:
- Find all the vertex points of the model
- Find the 'PointOn' of each beam segment
- Find the vector of each beam
- Determine whether the pointOn lies along the beam by using the cross product of the beam vector and the vector between the 'PointOn' point and the start point of the beam itself.
This is not only very slow, it does not work for some cases.
Does anyone have an alternative methodology for assigning beam section orientation to randomly oriented beams? Or, can anyone see what I have missed in my code to make it 100% robust?
I attach the code below - It is part of a much bigger script, but hopefully this is enough to go on.
Many Thanks
def orient():
import assembly
## Define some Abaqus shortcuts
a = mdb.models['Model-1'].rootAssembly
p = mdb.models['Model-1'].parts[final_part_name]
e = p.edges
##
def find_vectors():
## Find vectors of beams, calculated earlier in code
## and written to a file, make array of these for access now
global vectorstack
vectorstack = []
vectors_file = open('vectors.txt','r')
for lines in vectors_file:
lines = lines.rstrip('\n')
lines = separate(lines,',',0)
vectorstack.append(lines)
vectorstack = sorted(vectorstack)
vectors_file.close()
find_vectors()
print 'There are', len(vectorstack), 'vectors'
##
def find_internal_vertex_points():
## find number of vertices (all the points where the beams intersect)
## in model from Abaqus internal database
global list_of_vertex_points
list_of_vertex_points = []
vertex_no = len(mdb.models['Model-1'].parts['Merged_part'].vertices)
print 'There are', vertex_no, 'vertices'
## Append the vertex points to a list
for i in range(vertex_no):
## Loop through the list of vertices and extract each one
vertex = mdb.models['Model-1'].parts['Merged_part'].vertices[i]
## Remove the additional metadata attached to the coordinates
vertex_log = str(vertex)
vertex_log = vertex_log.lstrip(str("({'index': ")+str(i))
vertex_log = vertex_log.lstrip(str(", 'instanceName': None, 'pointOn':" ))
vertex_log = vertex_log.rstrip(',)})')
vertex_log = vertex_log.lstrip('((')
## Add the vertex point to list_of_vertex_points
vertex_log = vertex_log.split(',')
## Turn the points into an x,y,z list of floats
points = (float(vertex_log[0]),float(vertex_log[1]),float(vertex_log[2]))
list_of_vertex_points.append(points)
## Optionally, make a visible datum point at this point
## (useful for debugging)
#a.DatumPointByCoordinate(coords=(points))
## Sort the list of points
list_of_vertex_points = sorted(list_of_vertex_points)
print 'list_of_vertex_points', list_of_vertex_points
find_internal_vertex_points()
##
def find_reference_points_along_each_edge():
## Loop through the list of edges in Abaqus, this identifies
## the datum point associated with each edge, which can be used to
## select the edge and thus assign material properties
global the_ref_points_list
## Make blank list for coords of the edges
the_ref_points_list = []
## Determine how many edges there are
edges_no = len(mdb.models['Model-1'].parts['Merged_part'].edges)
print 'There are', edges_no, 'reference points along the edges'
##
for i in range(edges_no):
## Get data from Abaqus
this_edge = str(mdb.models['Model-1'].parts['Merged_part'].edges[i])
## Remove metadata and assign points as coords
this_edge = this_edge.split(',')
e3 = str(this_edge[len(this_edge)-2]).rstrip(')')
e2 = str(this_edge[len(this_edge)-3])
e1 = str(this_edge[len(this_edge)-4]).lstrip("'pointOn': ((")
## Determine the point along the edge
e1,e2,e3 = float(e1),float(e2),float(e3)
points = (e1,e2,e3)
## Add this to the list of edges
the_ref_points_list.append(points)
## Datum point along the edge
#a.DatumPointByCoordinate(coords=(points))
the_ref_points_list = sorted(the_ref_points_list)
find_reference_points_along_each_edge()
##
print 'list_of_edges: \n', the_ref_points_list
##
## Loop through all of the vertices
## Inside this, loop through all of the reference points
## Find the vector which connects the reference point to the vertex point [sub_beam]
## Inside this, loop through all of the beams, and see whether this sub_beam is parallel
## with an actual beam. If it is then actually asign the beam orientation
##
n = 0
list_of_assignments = []
## Looping thorugh all the vertex points
for vertices in list_of_vertex_points:
vx,vy,vz = vertices[0],vertices[1],vertices[2]
#a.DatumPointByCoordinate(coords=(vx,vy,vz))
#print vx,vy,vz
## Looping through all the reference points along each beam
for points in the_ref_points_list:
#print points
## Find the coordinates of the reference point
#a.DatumPointByCoordinate(coords=(points))
px,py,pz = points[0],points[1],points[2]
## find the vector between this reference point and each vertex
vecx,vecy,vecz = (px-vx),(py-vy),(pz-vz)
## find the vector of each beam
for beams in vectorstack:
bx,by,bz = float(beams[6]),float(beams[7]),float(beams[8])
#print 'bx,by,bz',bx,by,bz
## Determinants
deta = by*vecz - bz*vecy
detb = bz*vecx - bx*vecz
detc = bx*vecy - by*vecx
## Magnitude
cross_product = math.sqrt(deta**2+detb**2+detc**2)
#print 'cross_product', cross_product
## Moduli
mod1 = math.sqrt(bx**2 + by**2 + bz**2)
mod2 = math.sqrt(vecx**2 + vecy**2 + vecz**2)
#print 'Moduli', mod1, mod2
## Dot product
dot_product = (bx*vecx)+(by*vecy)+(bz*vecz)
#print 'Dot Product', dot_product
## Costheta
costheta = dot_product/(mod1*mod2)
#print 'Cos theta', costheta
## Determine whether to assign the properties by seeing whether
## the point vector and the beam vector are parallel
if cross_product == 0.0:
assign = True
elif costheta == 1.0:
assign = True
else:
assign = False
##
if assign == True:
n = n+1
# print 'beam is', bx,by,bz
# print 'sub beam is', vecx,vecy,vecz, '\n'
## Is beam 1,1,1, if so set random other vector
if abs(bx) == abs(by) and abs(bx) == abs(bz):
vec1,vec2,vec3 = 5.0,1.0,10.0
else:
vec1,vec2,vec3 = 1.0, 1.0,1.0
## Use vector above to calculate a normal
deta = by*vec3 - bz*vec2
detb = -(bz*vec1 - bx*vec3)
detc = bx*vec2 - by*vec1
## Assign section properties to the point
a.DatumPointByCoordinate(coords=(px,py,pz))
edges = e.findAt(((px,py,pz), ))
region=regionToolset.Region(edges=edges)
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(deta,-detb,detc))
orient()
[Non-text portions of this message have been removed]
I am trying to write a Python script to assign beam section orientation to a random network of beams. The part is made up of individual beams (in 3d space) and merged into a single part. Intersections are retained so that different parts can be assigned different properties. This causes an issue, as to assign beam section orientation to the part (necessary once it has been merged), it is necessary to select all of the regions between the intersections. If a single beam is intersected by one or more other beams, each part must be selected and an orientation assigned.
My current process for this is:
- Find all the vertex points of the model
- Find the 'PointOn' of each beam segment
- Find the vector of each beam
- Determine whether the pointOn lies along the beam by using the cross product of the beam vector and the vector between the 'PointOn' point and the start point of the beam itself.
This is not only very slow, it does not work for some cases.
Does anyone have an alternative methodology for assigning beam section orientation to randomly oriented beams? Or, can anyone see what I have missed in my code to make it 100% robust?
I attach the code below - It is part of a much bigger script, but hopefully this is enough to go on.
Many Thanks
def orient():
import assembly
## Define some Abaqus shortcuts
a = mdb.models['Model-1'].rootAssembly
p = mdb.models['Model-1'].parts[final_part_name]
e = p.edges
##
def find_vectors():
## Find vectors of beams, calculated earlier in code
## and written to a file, make array of these for access now
global vectorstack
vectorstack = []
vectors_file = open('vectors.txt','r')
for lines in vectors_file:
lines = lines.rstrip('\n')
lines = separate(lines,',',0)
vectorstack.append(lines)
vectorstack = sorted(vectorstack)
vectors_file.close()
find_vectors()
print 'There are', len(vectorstack), 'vectors'
##
def find_internal_vertex_points():
## find number of vertices (all the points where the beams intersect)
## in model from Abaqus internal database
global list_of_vertex_points
list_of_vertex_points = []
vertex_no = len(mdb.models['Model-1'].parts['Merged_part'].vertices)
print 'There are', vertex_no, 'vertices'
## Append the vertex points to a list
for i in range(vertex_no):
## Loop through the list of vertices and extract each one
vertex = mdb.models['Model-1'].parts['Merged_part'].vertices[i]
## Remove the additional metadata attached to the coordinates
vertex_log = str(vertex)
vertex_log = vertex_log.lstrip(str("({'index': ")+str(i))
vertex_log = vertex_log.lstrip(str(", 'instanceName': None, 'pointOn':" ))
vertex_log = vertex_log.rstrip(',)})')
vertex_log = vertex_log.lstrip('((')
## Add the vertex point to list_of_vertex_points
vertex_log = vertex_log.split(',')
## Turn the points into an x,y,z list of floats
points = (float(vertex_log[0]),float(vertex_log[1]),float(vertex_log[2]))
list_of_vertex_points.append(points)
## Optionally, make a visible datum point at this point
## (useful for debugging)
#a.DatumPointByCoordinate(coords=(points))
## Sort the list of points
list_of_vertex_points = sorted(list_of_vertex_points)
print 'list_of_vertex_points', list_of_vertex_points
find_internal_vertex_points()
##
def find_reference_points_along_each_edge():
## Loop through the list of edges in Abaqus, this identifies
## the datum point associated with each edge, which can be used to
## select the edge and thus assign material properties
global the_ref_points_list
## Make blank list for coords of the edges
the_ref_points_list = []
## Determine how many edges there are
edges_no = len(mdb.models['Model-1'].parts['Merged_part'].edges)
print 'There are', edges_no, 'reference points along the edges'
##
for i in range(edges_no):
## Get data from Abaqus
this_edge = str(mdb.models['Model-1'].parts['Merged_part'].edges[i])
## Remove metadata and assign points as coords
this_edge = this_edge.split(',')
e3 = str(this_edge[len(this_edge)-2]).rstrip(')')
e2 = str(this_edge[len(this_edge)-3])
e1 = str(this_edge[len(this_edge)-4]).lstrip("'pointOn': ((")
## Determine the point along the edge
e1,e2,e3 = float(e1),float(e2),float(e3)
points = (e1,e2,e3)
## Add this to the list of edges
the_ref_points_list.append(points)
## Datum point along the edge
#a.DatumPointByCoordinate(coords=(points))
the_ref_points_list = sorted(the_ref_points_list)
find_reference_points_along_each_edge()
##
print 'list_of_edges: \n', the_ref_points_list
##
## Loop through all of the vertices
## Inside this, loop through all of the reference points
## Find the vector which connects the reference point to the vertex point [sub_beam]
## Inside this, loop through all of the beams, and see whether this sub_beam is parallel
## with an actual beam. If it is then actually asign the beam orientation
##
n = 0
list_of_assignments = []
## Looping thorugh all the vertex points
for vertices in list_of_vertex_points:
vx,vy,vz = vertices[0],vertices[1],vertices[2]
#a.DatumPointByCoordinate(coords=(vx,vy,vz))
#print vx,vy,vz
## Looping through all the reference points along each beam
for points in the_ref_points_list:
#print points
## Find the coordinates of the reference point
#a.DatumPointByCoordinate(coords=(points))
px,py,pz = points[0],points[1],points[2]
## find the vector between this reference point and each vertex
vecx,vecy,vecz = (px-vx),(py-vy),(pz-vz)
## find the vector of each beam
for beams in vectorstack:
bx,by,bz = float(beams[6]),float(beams[7]),float(beams[8])
#print 'bx,by,bz',bx,by,bz
## Determinants
deta = by*vecz - bz*vecy
detb = bz*vecx - bx*vecz
detc = bx*vecy - by*vecx
## Magnitude
cross_product = math.sqrt(deta**2+detb**2+detc**2)
#print 'cross_product', cross_product
## Moduli
mod1 = math.sqrt(bx**2 + by**2 + bz**2)
mod2 = math.sqrt(vecx**2 + vecy**2 + vecz**2)
#print 'Moduli', mod1, mod2
## Dot product
dot_product = (bx*vecx)+(by*vecy)+(bz*vecz)
#print 'Dot Product', dot_product
## Costheta
costheta = dot_product/(mod1*mod2)
#print 'Cos theta', costheta
## Determine whether to assign the properties by seeing whether
## the point vector and the beam vector are parallel
if cross_product == 0.0:
assign = True
elif costheta == 1.0:
assign = True
else:
assign = False
##
if assign == True:
n = n+1
# print 'beam is', bx,by,bz
# print 'sub beam is', vecx,vecy,vecz, '\n'
## Is beam 1,1,1, if so set random other vector
if abs(bx) == abs(by) and abs(bx) == abs(bz):
vec1,vec2,vec3 = 5.0,1.0,10.0
else:
vec1,vec2,vec3 = 1.0, 1.0,1.0
## Use vector above to calculate a normal
deta = by*vec3 - bz*vec2
detb = -(bz*vec1 - bx*vec3)
detc = bx*vec2 - by*vec1
## Assign section properties to the point
a.DatumPointByCoordinate(coords=(px,py,pz))
edges = e.findAt(((px,py,pz), ))
region=regionToolset.Region(edges=edges)
p.assignBeamSectionOrientation(region=region, method=N1_COSINES, n1=(deta,-detb,detc))
orient()
[Non-text portions of this message have been removed]