Hello, I am trying to do the following via scripting:
- setup an Abaqus model, extracting the matrices
- computing in Python the eigenmodes
3. writing them to an odb
I have the following scripts, and I would like to know where the problem is, as I am getting different eigenmodes with respect to Abaqus jobs. Any input would be appreciated.
<i>from</i> abaqus <i>import</i> mdb, session
<i>from</i> abaqusConstants <i>import</i> *
<i>import</i> regionToolset
<i>import</i> mesh
<i>import</i> time
<i>import</i> numpy <i>as</i> np
model_name = 'BlockModel'
part_name = 'Block'
material_name = 'Steel'
section_name = 'BlockSection'
assembly_name = 'Assembly'
instance_name = part_name + '-1'
L = 1.0
H = 0.2
W = 0.2
seed_size = 0.1
E = 210e9
nu = 0.3
rho = 7850
cae_file = 'block_model.cae'
input_file = 'block_model.inp'
<i>if</i> model_name in mdb.models:
<i>del</i> mdb.models[model_name]
model = mdb.Model(<i>name</i>=model_name)
s = model.ConstrainedSketch(<i>name</i>='__profile__', <i>sheetSize</i>=2.0)
s.rectangle(<i>point1</i>=(0.0, 0.0), <i>point2</i>=(L, H))
part = model.Part(<i>name</i>=part_name, <i>dimensionality</i>=THREE_D, <i>type</i>=DEFORMABLE_BODY)
part.BaseSolidExtrude(<i>sketch</i>=s, <i>depth</i>=W)
material = model.Material(<i>name</i>=material_name)
material.Elastic(<i>table</i>=((E, nu),))
material.Density(<i>table</i>=((rho,),))
model.HomogeneousSolidSection(
<i>name</i>=section_name,
<i>material</i>=material_name,
<i>thickness</i>=None
)
cells = part.cells[:]
region = regionToolset.Region(<i>cells</i>=cells)
part.SectionAssignment(<i>region</i>=region, <i>sectionName</i>=section_name)
assembly = model.rootAssembly
instance = assembly.Instance(<i>name</i>=instance_name, <i>part</i>=part, <i>dependent</i>=ON)
faces = instance.faces
leftFaces = faces.getByBoundingBox(<i>xMin</i>=-1e-6, <i>xMax</i>=1e-6)
<i>if</i> len(leftFaces) == 0:
leftFaces = [f <i>for</i> f <i>in</i> faces <i>if</i> abs(f.pointOn[0]) < 1e-6]
assembly.Set(<i>name</i>='LeftFace', <i>faces</i>=leftFaces)
assembly.Set(<i>name</i>='WholeInstance', <i>nodes</i>=instance.nodes)
part.seedPart(<i>size</i>=seed_size, <i>deviationFactor</i>=0.1, <i>minSizeFactor</i>=0.1)
elemType1 = mesh.ElemType(<i>elemCode</i>=C3D8, <i>elemLibrary</i>=STANDARD)
part.setElementType(<i>regions</i>=(part.cells[:],), <i>elemTypes</i>=(elemType1,))
part.generateMesh()
assembly.regenerate()
leftSet = assembly.sets['LeftFace']
bc_labels = [node.label <i>for</i> node <i>in</i> leftSet.nodes]
np.savetxt('leftSet.txt', bc_labels)
<i># model.EncastreBC(name='BC_Clamped', createStepName='Initial', region=leftSet) # do not apply this BC here and lock it afterwards</i>
model.FrequencyStep(<i>name</i>='ExtractMatrices', <i>previous</i>='Initial', <i>numEigen</i>=1)
mdb.saveAs(<i>pathName</i>=cae_file)
mdb.Job(<i>name</i>='BlockMatrixJob', <i>model</i>=model_name).writeInput(<i>consistencyChecking</i>=OFF)
input_file = 'BlockMatrixJob.inp'
<i>with</i> open(input_file, 'r') <i>as</i> f:
lines = f.readlines()
target_line = "*Step, name=ExtractMatrices, nlgeom=NO, perturbation"
<i>for</i> i, line <i>in</i> enumerate(lines):
<i>if</i> target_line in line:
lines[i] = "*Step, name=ExtractMatrices\\n"
target_line = '*Frequency, eigensolver=Lanczos, sim, acoustic coupling=on, normalization=mass'
<i>for</i> i, line <i>in</i> enumerate(lines):
<i>if</i> target_line in line:
lines[i] = "*MATRIX GENERATE, STIFFNESS, MASS\\n"
target_line = "1, , , , ,"
<i>for</i> i, line <i>in</i> enumerate(lines):
<i>if</i> target_line in line:
lines[i] = "*MATRIX OUTPUT, STIFFNESS, MASS, FORMAT=MATRIX INPUT\\n"
<i>with</i> open(input_file, 'w') <i>as</i> f:
f.writelines(lines)
time.sleep(5)
mdb.JobFromInputFile(<i>name</i>='BlockMatrixJob', <i>inputFileName</i>=input_file)
mdb.jobs['BlockMatrixJob'].submit()
mdb.jobs['BlockMatrixJob'].waitForCompletion()<i>import</i> numpy <i>as</i> np
<i>import</i> matplotlib.pyplot <i>as</i> plt
<i>from</i> scipy.sparse <i>import</i> coo_matrix
<i>from</i> scipy.sparse.linalg <i>import</i> eigsh
bc_labels = np.loadtxt('leftSet.txt')
bc_labels = bc_labels.astype(int)-1
bc_dofs = []
n_dofs_per_node = 3
<i>for</i> label <i>in</i> bc_labels:
<i>for</i> dof <i>in</i> range(n_dofs_per_node):
bc_dofs.append(label * n_dofs_per_node + dof)
data = np.loadtxt('BlockMatrixJob_STIF1.mtx', <i>delimiter</i>=',')
node_i = data[:,0].astype(int) - 1
dof_i = data[:,1].astype(int) - 1
node_j = data[:,2].astype(int) - 1
dof_j = data[:,3].astype(int) - 1
vals = data[:,4]
n_dofs_per_node = 3
rows = node_i * n_dofs_per_node + dof_i
cols = node_j * n_dofs_per_node + dof_j
K = coo_matrix((vals, (rows, cols))).tocsr()
K = K.todense()
K = 1/2 * (K + K.T)
data = np.loadtxt('BlockMatrixJob_MASS1.mtx', <i>delimiter</i>=',')
node_i = data[:,0].astype(int) - 1
dof_i = data[:,1].astype(int) - 1
node_j = data[:,2].astype(int) - 1
dof_j = data[:,3].astype(int) - 1
vals = data[:,4]
n_dofs_per_node = 3
rows = node_i * n_dofs_per_node + dof_i
cols = node_j * n_dofs_per_node + dof_j
M = coo_matrix((vals, (rows, cols))).tocsr()
M = M.todense()
not_bc_dofs = np.setdiff1d(np.arange(K.shape[0]), bc_dofs)
Mr = M[np.ix_(not_bc_dofs, not_bc_dofs)]
Kr = K[np.ix_(not_bc_dofs, not_bc_dofs)]
vals, vecs = eigsh(Kr, 20, <i>M</i>=Mr)
freqs = np.sqrt(vals) / (2 * np.pi)
true_vecs = np.zeros((K.shape[0], vecs.shape[1]))
true_vecs[not_bc_dofs, :] = vecs
np.savetxt('eigenmodes.txt', true_vecs)
np.savetxt('eigenfreqs.txt', freqs)