...
 
Commits (2)
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import gdal
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import mesh_utilities as mesh
matplotlib.rcParams['xtick.direction']='out'
matplotlib.rcParams['ytick.direction']='out'
matplotlib.rcParams['axes.labelsize']='18'
matplotlib.rcParams['xtick.labelsize']='18'
matplotlib.rcParams['ytick.labelsize']='18'
matplotlib.rcParams['legend.fontsize']='18'
matplotlib.rcParams['font.family']='serif'
matplotlib.rcParams['font.sans-serif']='Times'
def get_map(natural_earth_file=None):
"""
zmap: elevation espessed in rgb scale 0-255 for plotting purposing only
usage: plt.imshow(zmap, extent=[-180,180,-90,90])
"""
if natural_earth_file == None:
data = "/work/tonini/data"
natural_earth_file = os.path.join(data,
'NE1_HR_LC_SR_W_DR/NE1_HR_LC_SR_W_DR.tif')
ds = gdal.Open(natural_earth_file)
data = ds.ReadAsArray()
zmap = np.moveaxis(data, 0, -1)
return zmap
def plot_mesh_abaqus(mesh_path, mesh_name):
mesh_d = mesh.read_abaqus(mesh_path, mesh_name)
polygons = mesh.create_polygons(mesh_d)
xmin = np.amin(polygons[:,:,0])
xmax = np.amax(polygons[:,:,0])
ymin = np.amin(polygons[:,:,1])
ymax = np.amax(polygons[:,:,1])
n_polygons = len(polygons[:,:,0])
plt.figure(figsize=(16,16))
plt.grid(True)
zmap = get_map()
plt.imshow(zmap, extent=[-180,180,-90,90])
for i in range(n_polygons):
plt.fill(np.append(polygons[i,:,0], polygons[i,0,0]),
np.append(polygons[i,:,1], polygons[i,0,1]),
fill=False, color="#333333", linewidth=1)
# plt.legend()
# plt.title(name)
plt.xlabel(r"Longitude ($^\circ$)")
plt.ylabel(r"Latitude ($^\circ$)")
plt.axis("equal")
plt.xlim(xmin-1, xmax+1)
plt.ylim(ymin-1, ymax+1)
fileout = os.path.join(mesh_path, mesh_name + "_mesh.png")
plt.savefig(fileout, fmt="png", bbox_inches="tight", dpi=150)
plt.close()
def main():
"""
"""
mesh_path = sys.argv[1]
mesh_name = sys.argv[2]
plot_mesh_abaqus(mesh_path, mesh_name)
if __name__ == "__main__":
main()
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import datetime
import os
import sys
import numpy as np
def abaqus2nodesfaces(mesh_path, mesh_name):
"""
"""
mesh_d = read_abaqus(mesh_path, mesh_name)
n_faces = len(mesh_d["faces"][:,0])
n_nodes = len(mesh_d["nodes"][:,0])
file_faces = os.path.join(mesh_path,
mesh_name + '_mesh_faces_test.dat')
with open(file_faces, "w") as f_faces:
for i in range(n_faces):
line = " ".join([str(item) for item in mesh_d["faces"][i,:]])
f_faces.write("{0}\n".format(line))
file_nodes = os.path.join(mesh_path,
mesh_name + '_mesh_nodes_test.dat')
with open(file_nodes, "w") as f_nodes:
for i in range(n_nodes):
ind = str(int(mesh_d["nodes"][i,0]))
line = " ".join([str(item) for item in mesh_d["nodes"][i,1:]])
f_nodes.write("{0} {1}\n".format(ind, line))
def nodesfaces2abaqus(mesh_path, mesh_name, name_out):
"""
"""
header_txt = ("*HEADING",
"this file is generated by pymesh.py: {0}".format(datetime.date.today()),
"version: no version",
"**",
"********************************** P A R T S **********************************",
"*PART, NAME=Part-Default",
"**",
"********************************** N O D E S **********************************",
"*NODE, NSET=ALLNODES")
middle_txt = ("**",
"********************************** E L E M E N T S ****************************",
"*ELEMENT, TYPE=STRI3, ELSET=EB1")
footer_txt = ("**",
"********************************** P R O P E R T I E S ************************",
"*SHELL SECTION, ELSET=EB1, SECTION INTEGRATION=SIMPSON, MATERIAL=Default-Steel",
"1.000000e+00",
"**",
"*END PART",
"**",
"**",
"**",
"********************************** E N D P A R T S **********************************",
"**",
"**",
"********************************** A S S E M B L Y ************************************",
"**",
"*ASSEMBLY, NAME=ASSEMBLY1",
"**",
"*INSTANCE, NAME=Part-Default_1, PART=Part-Default",
"*END INSTANCE",
"**",
"*END ASSEMBLY",
"**",
"**",
"**",
"*MATERIAL, NAME = Default-Steel",
"*ELASTIC, TYPE=ISOTROPIC",
"no value, no value",
"*DENSITY",
"no value",
"*CONDUCTIVITY,TYPE=ISO",
"no value",
"*SPECIFIC HEAT",
"no value",
"**",
"**",
"************************************** H I S T O R Y *************************************",
"**",
"*PREPRINT",
"**",
"**************************************** S T E P 1 ***************************************",
"*STEP,INC=100,NAME=Default Set",
"**",
"*STATIC",
"-, -, -, -",
"**",
"**",
"**",
"**",
"*END STEP")
mesh_txt = read_nodesfaces(mesh_path, mesh_name)
n_nodes = len(mesh_txt["nodes"])
n_elements = len(mesh_txt["faces"])
print(n_nodes, n_elements)
file_out = name_out + "_mesh.inp"
with open(file_out, "w") as f_out:
for h_line in header_txt:
f_out.write("{0}\n".format(h_line))
for i in range(n_nodes):
ind, lon, lat, z = mesh_txt["nodes"][i,:]
f_out.write("{:.0f}, {:f}, {:f}, {:f}\n".format(ind, lon, lat, z))
for m_line in middle_txt:
f_out.write("{0}\n".format(m_line))
for i in range(n_elements):
ind, i, j, k = mesh_txt["faces"][i,:]
f_out.write("{:.0f}, {:.0f}, {:.0f}, {:.0f}\n".format(ind, i, j, k))
for f_line in footer_txt:
f_out.write("{0}\n".format(f_line))
def ll2xyz():
"""
"""
pass
## utm zone (from epicentre)
# zone_n = hypocentre["utm_zone_number"]
# zone_l = hypocentre["utm_zone_letter"]
# mesh_d["nodes_utm"] = mesh_d["nodes"]
# n_nodes, tmp = np.shape(mesh_d["nodes"])
# tmp = ll2utm(mesh_d["nodes"][:, 1],
# mesh_d["nodes"][:, 2])
# mesh_d["nodes_utm"][:, 1] = tmp[0]
# mesh_d["nodes_utm"][:, 2] = tmp[1]
# for i in range(n_nodes):
# tmp = utm.from_latlon(mesh_d["nodes"][i,2],
# mesh_d["nodes"][i,1])
# mesh_d["nodes_utm"][i,1] = tmp[0]
# mesh_d["nodes_utm"][i,2] = tmp[1]
def read_nodesfaces(mesh_path, mesh_name):
"""
Read mesh file composed by the following text files:
- mesh_name_mesh_faces.dat
- mesh_name_mesh_nodes.dat
Return a python dictionary contaning name, faces and nodes
as keys of the dictionary.
"""
mesh_d = {}
file_mesh_faces = os.path.join(mesh_path,
mesh_name + '_mesh_faces.dat')
file_mesh_nodes = os.path.join(mesh_path,
mesh_name + '_mesh_nodes.dat')
mesh_d["name"] = mesh_name
mesh_d["faces"] = np.loadtxt(file_mesh_faces)
mesh_d["nodes"] = np.loadtxt(file_mesh_nodes)
return mesh_d
def read_abaqus(mesh_path, mesh_name, n_verts=None):
"""
"""
if n_verts is None:
n_verts = 3
mesh_d = {}
mesh_d["name"] = mesh_name
filename = os.path.join(mesh_path, mesh_name + "_mesh.inp")
nodes = []
with open(filename, "r") as f:
ic = 0
for line in f:
if (ic < 9):
pass
else:
if line.startswith("**"):
break
else:
tmp = line.strip().split(",")
inode = int(tmp[0])
nodes.append([inode] + [float(tmp[i]) for i in [1,2,3]])
ic += 1
newstart = ic + 3
faces = []
with open(filename, "r") as f:
ic = 0
for line in f:
if (ic < newstart):
pass
else:
if line.startswith("**"):
break
else:
tmp = line.strip().split(",")
faces.append([int(tmp[i]) for i in range(n_verts+1)])
ic += 1
mesh_d["nodes"] = np.array(nodes)
mesh_d["faces"] = np.array(faces)
return mesh_d
def create_polygons(mesh_dictionary):
"""
"""
faces = mesh_dictionary["faces"][:,1:].astype(int)
nodes = mesh_dictionary["nodes"][:,1:]
n_faces, n_vertices = np.shape(faces)
polygons = np.zeros((n_faces, n_vertices, 3))
for ip in range(n_faces):
ind = faces[ip,:]
for iv in range(n_vertices):
polygons[ip,iv,:] = nodes[ind[iv]-1,:]
return polygons
def main():
"""
The script can be run as a standalone tool.
Example:
python abaqus.py /scratch/projects/cat/SPTHA/EXERCISES/MAKRAN/data/ makran
"""
mesh_path = sys.argv[1]
mesh_name = sys.argv[2]
m = read_nodesfaces(mesh_path, mesh_name)
#file_out = os.path.join(mesh_path, mesh_name)
#mesh_abaqus_tmp = nodesfaces2abaqus(mesh_path, mesh_name, file_out)
abaqus2nodesfaces(mesh_path, mesh_name)
print(np.shape(create_polygons(m)))
if __name__ == "__main__":
main()