"""Analyze the results from a wheel substructure and create the
necessary data structures to setup the user element.
.. codeauthor:: Knut Andreas Meyer
"""
# Python imports
from __future__ import print_function
import numpy as np
# Abaqus imports
from abaqusConstants import *
# Project imports
import rollover.utils.abaqus_python_tools as apt
import rollover.utils.naming_mod as names
[docs]def get_uel_mesh(quadratic_elements=True):
"""Determine the mesh from the substructure simulation.
Produces the following files:
- `names.uel_stiffness_file`: The stiffness matrix to be read by the
fortran uel subroutine
- `names.uel_coordinates_file`: The coordinates of the contact nodes
in the user element.
- `names.uel_elements_file`: The indices of the user element nodes
that belong to each element.
"""
ke_raw = get_stiffness(names.substr_mtx_file)
rp_nr, contact_node_labels = get_mtx_nodes(names.substr_mtx_file)
ke = reorder_stiffness(ke_raw, rp_nr)
coords = get_node_coords(names.substr_node_coords_file,
names.substr_node_labels_file,
contact_node_labels)
if quadratic_elements:
elements = get_element_connectivity_quad(coords)
else:
elements = get_element_connectivity(coords)
save_uel(ke, coords, elements)
[docs]def get_stiffness(mtx_file):
"""Extracts the stiffness from the mtx file `mtx_file`, which was
generated by an Abaqus substructure generate step.
:param mtx_file: Name of the mtx file to read the stiffness matrix
from.
:type mtx_file: str
:returns: The stiffness matrix
:rtype: np.array
"""
with open(mtx_file, 'r') as mtx:
mtx_str = mtx.read()
mat_str = mtx_str.split('*MATRIX,TYPE=STIFFNESS')[-1].split('*')[0].strip(',').strip('\n')
mat_vec = []
for entry in mat_str.split():
ent = entry.strip(',').strip('\n')
try:
mat_vec.append(float(ent))
except ValueError as e:
if len(ent) == 0:
pass
else:
print('Cannot convert "' + ent + '" to a float')
raise e
mat_vec = np.array(mat_vec)
ndof = -0.5+np.sqrt(0.25+mat_vec.size*2)
if np.abs(ndof-int(ndof)) < 1.e-10:
ndof = int(ndof)
else:
print('Error reading matrix from ' + mtx_file + '.mtx')
return None
kmat = np.zeros((ndof,ndof))
k = 0
for i in range(ndof):
for j in range(i+1):
kmat[i,j] = mat_vec[k]
kmat[j,i] = kmat[i,j]
k = k + 1
return kmat
[docs]def get_mtx_nodes(mtx_file):
"""Extracts the node labels from the mtx file `mtx_file`, which was
generated by an Abaqus substructure generate step. Note, node
numbers starts from zero, while node labels starts from 1.
:param mtx_file: Name of the mtx file to read the node labels matrix
from.
:type mtx_file: str
:returns: List with items
- Element node index for reference point node (int)
(Has to be between 0 and number of nodes in element)
- np.array of contact node labels corresponding to the
node labels in names.uel_contact_node_labels_file.
(These are not restricted to be less than the number of
nodes in the element)
:rtype: list
"""
with open(mtx_file, 'r') as mtx:
# Skip the introduction
line = mtx.readline()
while not line.startswith('** ELEMENT NODES'):
line = mtx.readline()
# Read the element nodes category
line = mtx.readline()
node_str = ''
while line.startswith('**'):
node_str = node_str + line[2:]
line = mtx.readline()
node_inds = [int(n) for n in node_str.split(',')]
# Read the node dofs to determine which node is the reference
# point
node_dofs = []
node_dofs_ind = 0
while not line.startswith('*'):
node_dofs.append([int(s) for s in line.split(',')])
if len(node_dofs[-1]) == 7: # All disp and rotations
rp_dof_ind = node_dofs_ind
node_dofs_ind += 1
line = mtx.readline()
rp_node_nr = node_dofs[rp_dof_ind][0] - 1
contact_nodes = [node_ind for i, node_ind in enumerate(node_inds) if i != rp_node_nr]
return rp_node_nr, contact_nodes
[docs]def reorder_stiffness(ke_raw, rp_nr):
"""Reorder the stiffness such that the dofs related to the reference
point comes first. The order of the remaining nodes is unaffacted.
:param ke_raw: Unordered stiffness matrix
:type ke_raw: np.array
:param rp_nr: Node number in the element for the reference point
:type rp_nr: int
:returns: Ordered stiffness matrix
:rtype: np.array
"""
ndof_trans = 3
ndof_rot = 3
ndof = ke_raw.shape[0]
nnods = (ndof - ndof_rot)/ndof_trans # Number of nodes incl. rp
# Add the dofs for the reference point first
reorder = [ndof_trans*rp_nr + i for i in range(ndof_trans+ndof_rot)]
# Then add all nodes that previously were before the rp
for node_nr in range(rp_nr):
for dof_nr in range(ndof_trans):
reorder.append(node_nr*ndof_trans + dof_nr)
# Finally add nodes that previously were after the rp
for node_nr in range(rp_nr+1, nnods):
for dof_nr in range(ndof_trans):
reorder.append(ndof_rot + node_nr*ndof_trans + dof_nr)
reorder = np.array(reorder, dtype=np.int)
return ke_raw[np.ix_(reorder, reorder)]
[docs]def get_node_coords(coords_file, labels_file=None, contact_node_labels=None):
""" {TEST}
Get coordinates from `coords_file` with labels according to
`labels_file`. The coordinates are sorted such that the labels match
`contact_node_labels` if this list is present.
:param coords_file: .npy file containing node coordinates
:type coords_file: str
:param labels_file: .npy file containing node labels. Required if
`contact_node_labels` is not None.
:type labels_file: str
:param contact_node_labels: List of contact node labels that may
have a different order than those from
the `labels_file`. The returned
coordinate list will be sorted to follow
the order of contact_node_labels, if not
None. Defaults to None.
:type contact_node_labels: iterable[ int ]
:returns: The coordinates from `coords_file`, sorted according to
`contact_node_labels`, if it is present.
:rtype: np.array
"""
coords = np.load(coords_file)
if contact_node_labels is not None:
labels = list(np.load(labels_file))
sort_inds = np.array([labels.index(label) for label in contact_node_labels], dtype=np.int)
coords = coords[sort_inds]
return coords
[docs]def get_element_connectivity(coords):
""" Knowing that the mesh is revolved around the x-axis and that
we only have coordinates of the contact nodes. Create the element
connectivity (which nodes belong to which element) for the quoad
mesh.
:param coords: List of x,y,z coordinates generated by revolution
around the x-axis
:type coords: np.array (shape = [npoints, 3])
:returns: List of which 4 node indices that belong to each element.
:rtype: list[ list[ int ] ]
"""
# Create an index matrix where the row goes along the section and
# the column along the angular direction.
index_matrix = get_mesh_inds(coords)
elems = []
for row1, row2 in zip(index_matrix[:-1], index_matrix[1:]):
for inds in zip(row1[:-1], row1[1:], row2[1:], row2[:-1]):
elems.append(inds)
return elems
[docs]def get_element_connectivity_quad(coords):
""" Knowing that the mesh is revolved around the x-axis and that
we only have coordinates of the contact nodes. Create the element
connectivity (which nodes belong to which element) for the quoad
mesh.
:param coords: List of x,y,z coordinates generated by revolution
around the x-axis
:type coords: np.array (shape = [npoints, 3])
:returns: List of which 4 node indices that belong to each element.
:rtype: list[ list[ int ] ]
"""
# Create an index matrix where the row goes along the section and
# the column along the angular direction.
index_matrix = get_mesh_inds(coords)
elems = []
for row1, row2, row3 in zip(index_matrix[:-2:2], index_matrix[1:-1:2], index_matrix[2::2]):
for inds in zip(row1[:-2:2], row1[2::2], row3[2::2], row3[:-2:2],
row1[1:-1:2], row2[1::1], row3[1:-1:2], row2[:-1:1]):
# Element numbering for 8 node membrane element
# row1: 1 5 2
# row2: 8 6
# row3: 4 7 3
elems.append(inds)
return elems
[docs]def get_mesh_inds(coords):
""" {TEST}
Given a list of randomly sorted coordinates `coords` representing
points generated by a revolution pattern from points on a curve in
the xy-plane around the x-axis: Determine an index matrix such the
indices of the nearest neighbours can be located via the matrix.
:Note: Limitations of the current implementation:
- The points on the initial curve must be sufficiently spaced
in the x-direction.
- A full revolution is not supported and only points crossing
the xy-plane with a negative y-coordinate is supported.
:param coords: List of x,y,z coordinates generated by revolution
around the x-axis
:type coords: np.array (shape = [npoints, 3])
:returns: A matrix with indices corresponding to the elements of
`coords`. The first index goes in the positive angular
direction around the x-axis, while the second goes in the
positive x-direction.
:rtype: list[ list[ int ] ]
"""
TOL = 1.e-2 # Linear tolerance (length unit)
# Get angle around the x-axis, measured from the negative y-axis
angles = np.arctan2(-coords[:, 2], -coords[:, 1])
# Get radius and x-position
radii = np.sqrt(coords[:, 1]**2 + coords[:, 2]**2)
xcoords = coords[:, 0]
ang_tol = TOL/np.max(radii)
unique_angles = get_unique(angles, ang_tol)
unique_xcoords = get_unique(xcoords, TOL)
index_matrix = []
for ang in unique_angles:
index_matrix.append([])
last_failed = True # The first x-coordinate should be found
for xcoord in unique_xcoords:
try:
coord_index = find_coord(find_coords=(ang, xcoord),
search_coords=(angles, xcoords),
tol=[ang_tol, TOL])
this_failed = False
except ValueError:
this_failed = True
if last_failed and this_failed:
raise ValueError('Could not determine coordinates')
elif not this_failed:
index_matrix[-1].append(coord_index)
last_failed = this_failed
return index_matrix
[docs]def get_unique(vector, tol=0):
""" {TEST}
Given a vector `vector`, return a new vector containing only the
unique entries in vector. The output is sorted. If a tolerance is
given, it represents the maximum difference between two elements
considered to be non-unique.
:param vector: Iterable from which unique values will be identified.
:type vector: iterable[ float / int ]
:param tol: Tolerance within which elements of `vector` should be
considered unique, defaults to 0
:type tol: float / int
:returns: List of unique (within tol) values in `vector`
:rtype: np.array
"""
sorted = np.sort(vector)
unique_vals = []
current_vals = [sorted[0]]
for v in sorted[1:]:
if v > current_vals[0] + tol:
unique_vals.append(np.average(current_vals))
current_vals = [v]
else:
current_vals.append(v)
unique_vals.append(np.average(current_vals))
return np.array(unique_vals)
[docs]def find_coord(find_coords, search_coords, tol=1.e-6):
""" {TEST}
Find the index for the coordinate in search_coords that matches
the coordinate find_coords.
:param find_coords: Coordinates for the point to find
:type find_coords: tuple[ float ]
:param search_coords: Coordinate lists to be searched through for
match. Length of tuple must match
`find_coords`
:type search_coords: tuple[ np.array ]
:param tol: Tolerance for the found coordinate to be considered a
match. If tuple, the length must match `find_coords`
:type tol: float / tuple[ float ]
:returns: Index of the found coordinate
:rtype: int
"""
if isinstance(tol, float):
tol_list = [tol for _ in find_coords]
else:
tol_list = tol[:]
dist2 = 0.0
for find_coord, search_coord, the_tol in zip(find_coords, search_coords, tol_list):
dist2 = dist2 + ((find_coord - search_coord)/the_tol)**2
min_ind = np.argmin(dist2)
if dist2[min_ind] > 1.0:
raise ValueError('Could not identify matching coordinate within tol, '
+ 'error is %0.2f %% of tol' % (100.0*np.sqrt(dist2[min_ind])))
return min_ind
[docs]def save_uel(stiffness, coordinates, elements):
""" Save the stiffness, node coordinates and element connectivity
for the user element to be imported. Stiffness will be read by
fortran subroutine, while coordinates and elements will be read by
abaqus python when setting up the new simulation.
:param stiffness: Stiffness matrix, will be saved to
`names.uel_stiffness_file`
:type stiffness: np.array
:param coordinates: Node coordinates, will be saved to
`names.uel_coordinates_file`
:type coordinates: np.array
:param elements: Element connectivity (nodes in coordinates
belonging to which element), will be saved to
`names.uel_elements_file`
:type elements: np.array
:returns: None
:rtype: None
"""
# Create file to import stiffness matrix in fortran uel subroutine
with open(names.uel_stiffness_file, 'w') as fid:
ndof = stiffness.shape[0]
fid.write('%5u\n' % ndof) # First line for allocating matrix
for i in range(ndof):
for j in range(i, ndof):
fid.write('%25.15e\n' % stiffness[i,j])
# If verifying that correct indices are transferred,
# use the following format:
# fid.write('%5u, %5u, %25.15e\n' % (i+1, j+1, stiffness[i,j]))
# Create file with node coordinates
np.save(file=names.uel_coordinates_file, arr=coordinates)
# Create file with element nodes
np.save(file=names.uel_elements_file, arr=elements)
[docs]def create_test_part(quadratic_elements=True):
""" Create a test part to verify that the elements and nodes are
identified correctly
"""
if quadratic_elements:
elem_shape=QUAD8
else:
elem_shape=QUAD4
coords = np.load('uel_coordinates.npy')
element_connectivity = np.load('uel_elements.npy')
wheel_model = apt.create_model('WHEEL_TEST_IMP')
wheel_part = wheel_model.Part(name=names.wheel_part, dimensionality=THREE_D,
type=DEFORMABLE_BODY)
nodes = [wheel_part.Node(coordinates=coord) for coord in coords]
elems = []
for ec in element_connectivity:
enodes = [nodes[i] for i in ec]
elems.append(wheel_part.Element(nodes=enodes, elemShape=elem_shape))