plots.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#####################################################################################
# This file is part of NS2DDV. #
# #
# Copyright(C) 2011-2018 C. Calgaro (caterina.calgaro@math.univ-lille1.fr) #
# E. Creusé (emmanuel.creuse@math.univ-lille1.fr) #
# T. Goudon (thierry.goudon@inria.fr) #
# A. Mouton (alexandre.mouton@math.univ-lille1.fr) #
# #
# NS2DDV is free software: you can redistribute it and/or modify it under the terms #
# of the GNU General Public License as published by the Free Software Foundation, #
# either version 3 of the License, or (at your option) any later version. #
# #
# NS2DDV is distributed in the hope that it will be useful, but WITHOUT ANY #
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A #
# PARTICULAR PURPOSE. See the GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License along with #
# NS2DDV. If not, see <http://www.gnu.org/licenses/>. #
#####################################################################################
from base import *
#from readlog import *
#import psutil
def plot_mesh(meshh5, method, aliasfield='', \
list_bbox=['default']*4, \
tuple_inches_figsize=(12.,12.), \
outputpng='./temp.png', \
forcesave=False):
if (aliasfield == ''):
newaliasfield = 'Mesh'
else:
newaliasfield = aliasfield
triangulation = build_triangulation(meshh5, method)
f = h5py.File(meshh5, 'r')
geom_kind = f.attrs['Geometry'].decode('UTF-8')
f.close()
fs = 20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
matplotlib.rcParams['font.family'] = 'monospace'
fig = plt.figure(figsize=tuple_inches_figsize)
ax = plt.gca()
ax.set_aspect('equal')
im = plt.triplot(triangulation)
xnodes = triangulation.x
ynodes = triangulation.y
xmin = min(xnodes)
xmax = max(xnodes)
ymin = min(ynodes)
ymax = max(ynodes)
if (geom_kind == 'DIHEDRON'):
triangles = triangulation.triangles
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yA1_index = i
themin = ynodes[i]
yA1 = ynodes[yA1_index]
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xA2_index = i
themin = xnodes[i]
xA2 = xnodes[xA2_index]
xslopemin = xA2
yslopemin = ymax
for i in range(len(triangles)):
if (xA2_index in triangles[i]):
x1 = xnodes[triangles[i][0]]
y1 = ynodes[triangles[i][0]]
x2 = xnodes[triangles[i][1]]
y2 = ynodes[triangles[i][1]]
x3 = xnodes[triangles[i][2]]
y3 = ynodes[triangles[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles[i][2]
xA1 = (yA1-ymin)*(xslopemin-xA2)/(yslopemin-ymin) + xA2
xs = [xmin, xA1, xA2, xmin]
ys = [yA1, yA1, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
elif (geom_kind == 'STEP'):
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminleft_index = i
themin = ynodes[i]
yminleft = ynodes[yminleft_index]
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmax)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminright_index = i
themin = ynodes[i]
yminright = ynodes[yminright_index]
if (abs(yminleft-ymin)/(ymax-ymin) < 1e-8):
ystep = yminright
# Step is on the right part of the domain
themax = xmin
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] >= themax):
xstep_index = i
themax = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xstep, xmax, xmax, xstep]
ys = [ystep, ystep, ymin, ymin]
elif (abs(yminright-ymin)/(ymax-ymin) < 1e-8):
ystep = yminleft
# Step is on the left part of the domain
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xstep_index = i
themin = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xmin, xstep, xstep, xmin]
ys = [ystep, ystep, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
plt.xlabel("x", fontsize=fs)
plt.ylabel("y", fontsize=fs)
plt.title(newaliasfield, fontsize=fs)
ttl = ax.title
ttl.set_position([.5, 1.1])
if (len(list_bbox) == 4):
if (list_bbox[0] == 'default'):
xminset = xmin
else:
xminset = list_bbox[0]
if (list_bbox[1] == 'default'):
xmaxset = xmax
else:
xmaxset = list_bbox[1]
if (list_bbox[2] == 'default'):
yminset = ymin
else:
yminset = list_bbox[2]
if (list_bbox[3] == 'default'):
ymaxset = ymax
else:
ymaxset = list_bbox[3]
ax.set_xlim([xminset, xmaxset])
ax.set_ylim([yminset, ymaxset])
else:
print('Error: list_bbox must be a list of size 4')
exit()
plt.draw()
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with file " + outputpng + ": the file already exists. Please choose another name")
exit()
if (not(outputpng == '')):
fig.savefig(outputpng, bbox_extra_artists=(ttl,), bbox_inches='tight')
del triangulation
if (geom_kind == 'DIHEDRON'):
del xnodes, ynodes, triangles#, x, y1, y2
fig.clf()
plt.close()
def plot_contour_withtri(inputh5, field, aliasfield, levels, list_bbox, triangulation, \
tuple_anchor_legend, tuple_inches_figsize, outputpng):
h5file = h5py.File(inputh5, 'r')
geom_kind = h5file.attrs['Geometry'].decode('UTF-8')
if (field in h5file.keys()):
data = h5file[field].value
else:
print('Error: cannot find the field '+field+' in file '+inputh5)
exit()
data = h5file[field].value
if ((not(data.shape[0] == 1)) and (data.shape[1] == 1)):
# Need to transpose the data
data = numpy.reshape(data, (1,data.shape[0]))[0]
else:
data = data[0]
min_data = min(data)
max_data = max(data)
#print([min_data,max_data])
#print([min(levels), max(levels)])
if ('Time' in h5file.keys()):
time = h5file['Time'].value[0]
else:
print('Error: cannot find the field ''Time'' in file '+inputh5)
exit()
if ('Reynolds_number' in h5file.keys()):
Re = h5file['Reynolds_number'].value[0]
else:
print('Error: cannot find the field ''Reynolds_number'' in file '+inputh5)
exit()
if ((max_data-min_data) > 1e-12*max([abs(max_data),abs(min_data)])):
constant = False
newlevels = levels.copy()
newlevels[0] += (max(levels)-min(levels))*0.000001
newlevels[-1] -= (max(levels)-min(levels))*0.000001
labels = []
for i in range(len(newlevels)):
if (newlevels[-1-i] < 0):
labels.append(aliasfield+" = {:.5e}".format(newlevels[-1-i]))
else:
labels.append(aliasfield+" = {:.5e} ".format(newlevels[-1-i]))
else:
constant = True
labels = [aliasfield+" = {:.5e}".format(min_data), "(Constant)"]
newlevels = []
fs = 20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
matplotlib.rcParams['font.family'] = 'monospace'
fig = plt.figure(figsize=tuple_inches_figsize)
ax = plt.gca()
ax.set_aspect('equal')
if (constant):
if (len(data) == len(triangulation.x)):
im = ax.tricontourf(triangulation.x, triangulation.y, data, colors='w')
elif (len(data) == len(triangulation.xtricenter)):
im = ax.tricontourf(triangulation.xtricenter, triangulation.ytricenter, data, colors='w')
for i in range(len(labels)):
im.collections[i].set_label(labels[i])
im.ax.tick_params(direction='out')
else:
if (len(data) == len(triangulation.x)):
im = ax.tricontour(triangulation.x, triangulation.y, data, levels=newlevels, extend='both', \
linewidths=3, cmap=plt.cm.jet)
elif (len(data) == len(triangulation.xtricenter)):
im = ax.tricontour(triangulation.xtricenter, triangulation.ytricenter, data, levels=newlevels, \
extend='both', \
linewidths=3, cmap=plt.cm.jet)
im.set_clim(min(levels),max(levels))
im.ax.tick_params(direction='out')
col = matplotlib.cm.get_cmap('jet')
codecolors = []
newcollection = []
for i in range(len(newlevels)):
thelevelrenorm = (newlevels[i]-min(newlevels))/(max(newlevels)-min(newlevels))
codecolor = col(thelevelrenorm)
codecolors.append(list(codecolor[0:3]))
thepatch = matplotlib.lines.Line2D([], [], color=codecolor, linestyle='solid', linewidth=3, label=labels[i])
newcollection.append(thepatch)
#a = dir(im.collections[0])
#for i in range(len(a)):
# print(a[i])
#print(im.collections[0].get_color())
#print(im.collections[0].get_offsets())
#exit()
#print(im.collections[0].get_array())
#print(len(im.collections))
ll = im.get_array()
for i in range(0,len(im.collections)):
for j in range(0,len(levels)):
if (abs(ll[i]-newlevels[j]) < 1e-12*abs(max(newlevels)-min(newlevels))):
im.collections[i].set_label(labels[j])
break
xnodes = triangulation.x
ynodes = triangulation.y
xmin = min(xnodes)
xmax = max(xnodes)
ymin = min(ynodes)
ymax = max(ynodes)
if (geom_kind == 'DIHEDRON'):
triangles = triangulation.triangles
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yA1_index = i
themin = ynodes[i]
yA1 = ynodes[yA1_index]
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xA2_index = i
themin = xnodes[i]
xA2 = xnodes[xA2_index]
xslopemin = xA2
yslopemin = ymax
for i in range(len(triangles)):
if (xA2_index in triangles[i]):
x1 = xnodes[triangles[i][0]]
y1 = ynodes[triangles[i][0]]
x2 = xnodes[triangles[i][1]]
y2 = ynodes[triangles[i][1]]
x3 = xnodes[triangles[i][2]]
y3 = ynodes[triangles[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles[i][2]
xA1 = (yA1-ymin)*(xslopemin-xA2)/(yslopemin-ymin) + xA2
xs = [xmin, xA1, xA2, xmin]
ys = [yA1, yA1, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
elif (geom_kind == 'STEP'):
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminleft_index = i
themin = ynodes[i]
yminleft = ynodes[yminleft_index]
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmax)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminright_index = i
themin = ynodes[i]
yminright = ynodes[yminright_index]
if (abs(yminleft-ymin)/(ymax-ymin) < 1e-8):
ystep = yminright
# Step is on the right part of the domain
themax = xmin
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] >= themax):
xstep_index = i
themax = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xstep, xmax, xmax, xstep]
ys = [ystep, ystep, ymin, ymin]
elif (abs(yminright-ymin)/(ymax-ymin) < 1e-8):
ystep = yminleft
# Step is on the left part of the domain
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xstep_index = i
themin = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xmin, xstep, xstep, xmin]
ys = [ystep, ystep, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
plt.xlabel("x", fontsize=fs)
plt.ylabel("y", fontsize=fs)
plt.title(aliasfield+" at time t = {:f}".format(time)+" (Re = {:f})".format(Re), fontsize=fs)
ttl = ax.title
ttl.set_position([.5, 1.1])
if (len(list_bbox) == 4):
if (list_bbox[0] == 'default'):
xminset = xmin
else:
xminset = list_bbox[0]
if (list_bbox[1] == 'default'):
xmaxset = xmax
else:
xmaxset = list_bbox[1]
if (list_bbox[2] == 'default'):
yminset = ymin
else:
yminset = list_bbox[2]
if (list_bbox[3] == 'default'):
ymaxset = ymax
else:
ymaxset = list_bbox[3]
ax.set_xlim([xminset, xmaxset])
ax.set_ylim([yminset, ymaxset])
else:
print('Error: list_bbox must be a list of size 4')
exit()
[oldhandles, labels] = im.ax.get_legend_handles_labels()
if (not(constant)):
oldhandles = newcollection
separator_patch = matplotlib.lines.Line2D([], [], color='w', label='')
leg_patch = matplotlib.lines.Line2D([], [], color='black', linestyle='solid', linewidth=3, label=aliasfield)
if (min_data < 0):
min_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e}'.format(min_data))
else:
min_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e} '.format(min_data))
if (max_data < 0):
max_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e}'.format(max_data))
else:
max_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e} '.format(max_data))
oldhandles.append(separator_patch)
oldhandles.append(leg_patch)
oldhandles.append(max_patch)
oldhandles.append(min_patch)
lgd = plt.legend(handles=oldhandles, fontsize=fs, ncol=1, loc='upper left', \
frameon=False, bbox_to_anchor=tuple_anchor_legend)
plt.draw()
if (not(outputpng == '')):
fig.savefig(outputpng, bbox_extra_artists=(lgd,ttl,), bbox_inches='tight')
del data, newlevels, xnodes, ynodes
if (geom_kind == 'DIHEDRON'):
del triangles#, x, y1, y2
fig.clf()
plt.close()
h5file.close()
def plot_contour(inputh5, field, \
aliasfield='', \
levels=[], \
list_bbox=['default']*4, \
tuple_anchor_legend=(1.01,1.01), \
tuple_inches_figsize=(12., 12.), \
outputpng='./temp.png', \
forcesave=False):
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
f = h5py.File(inputh5, 'r')
meshh5 = f.attrs['Mesh'].decode('UTF-8')
if (field in list_Pp):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif (field in list_Pv):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif (field in list_Prho):
method = f.attrs['Method_density'].decode('UTF-8')
f.close()
if (len(levels) == 0):
# Need to provide isovalues
theminmax = extract_minmax(inputh5, field)
newlevels = numpy.linspace(min(theminmax), max(theminmax), 10, endpoint=True)
else:
newlevels = levels.copy()
if (aliasfield == ''):
newaliasfield = field
else:
newaliasfield = aliasfield
# We verify if the mesh file is correct
if (not(os.path.isfile(meshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + meshh5.split('/')[-1]
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
meshh5 = testmeshh5
# Load mesh
triangulation = build_triangulation(meshh5, method)
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with file " + outputpng + ": the file already exists. Please choose another name")
exit()
plot_contour_withtri(inputh5, field, newaliasfield, newlevels, list_bbox, triangulation, \
tuple_anchor_legend, tuple_inches_figsize, outputpng)
def plot_two_contour_withtri(inputh5, field, aliasfield, list_alias_diff, levels, list_bbox, \
triangulations, tuple_anchor_legend, tuple_inches_figsize, outputpng):
h5file_1 = h5py.File(inputh5[0], 'r')
h5file_2 = h5py.File(inputh5[1], 'r')
geom_kind_1 = h5file_1.attrs['Geometry'].decode('UTF-8')
geom_kind_2 = h5file_2.attrs['Geometry'].decode('UTF-8')
if (geom_kind_1 == geom_kind_2):
geom_kind = geom_kind_1
else:
print('Error: the input HDF5 are not based on the same domain geometry kind')
exit()
if (field in h5file_1.keys()):
data_1 = h5file_1[field].value
else:
print('Error: cannot find the field '+field+' in file '+inputh5[0])
exit()
if ((not(data_1.shape[0] == 1)) and (data_1.shape[1] == 1)):
# Need to transpose the data
data_1 = numpy.reshape(data_1, (1,data_1.shape[0]))[0]
else:
data_1 = data_1[0]
if (field in h5file_2.keys()):
data_2 = h5file_2[field].value
else:
print('Error: cannot find the field '+field+' in file '+inputh5[1])
exit()
if ((not(data_2.shape[0] == 1)) and (data_2.shape[1] == 1)):
# Need to transpose the data
data_2 = numpy.reshape(data_2, (1,data_2.shape[0]))[0]
else:
data_2 = data_2[0]
min_data_1 = min(data_1)
max_data_1 = max(data_1)
min_data_2 = min(data_2)
max_data_2 = max(data_2)
if ('Time' in h5file_1.keys()):
time_1 = h5file_1['Time'].value[0]
else:
print('Error: cannot find the field ''Time'' in file '+inputh5[0])
exit()
if ('Time' in h5file_2.keys()):
time_2 = h5file_2['Time'].value[0]
else:
print('Error: cannot find the field ''Time'' in file '+inputh5[1])
exit()
if (abs(time_1-time_2) > 1e-12*max([time_1,time_2])):
print('Error: the following files do not contains data at same time step')
print(' '+inputh5[0]+' contains data at time t = '+str(time_1))
print(' '+inputh5[1]+' contains data at time t = '+str(time_2))
exit()
else:
time = time_1
if ('Reynolds_number' in h5file_1.keys()):
Re_1 = h5file_1['Reynolds_number'].value[0]
else:
print('Error: cannot find the field ''Reynolds_number'' in file '+inputh5[0])
exit()
if ('Reynolds_number' in h5file_2.keys()):
Re_2 = h5file_2['Reynolds_number'].value[0]
else:
print('Error: cannot find the field ''Reynolds_number'' in file '+inputh5[1])
exit()
if (abs(Re_1-Re_2) > 1e-12*max([Re_1,Re_2])):
print('Error: the following files do not contains data at same Reynolds number')
print(' '+inputh5[0]+' contains data with Re = '+str(Re_1))
print(' '+inputh5[1]+' contains data with Re = '+str(Re_2))
exit()
else:
Re = Re_1
if (((max_data_1-min_data_1) > 1e-12*max([abs(max_data_1),abs(min_data_1)])) \
and ((max_data_2-min_data_2) > 1e-12*max([abs(max_data_2),abs(min_data_2)]))):
constant = False
newlevels = levels.copy()
newlevels[0] += (newlevels[-1]-newlevels[0])*0.000001
newlevels[-1] -= (newlevels[-1]-newlevels[0])*0.000001
labels = []
for i in range(len(newlevels)):
if (newlevels[-1-i] < 0):
labels.append(aliasfield+" = {:.5e}".format(newlevels[-1-i]))
else:
labels.append(aliasfield+" = {:.5e} ".format(newlevels[-1-i]))
else:
constant = True
labels = [aliasfield+" = {:.5e}".format(min_data_1), "("+list_alias_diff[0]+" - Constant)", \
aliasfield+" = {:.5e}".format(min_data_2), "("+list_alias_diff[1]+" - Constant)"]
newlevels = []
fs = 20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
matplotlib.rcParams['font.family'] = 'monospace'
fig = plt.figure(figsize=tuple_inches_figsize)
ax = plt.gca()
ax.set_aspect('equal')
if (constant):
if (len(data_2) == len(triangulations[1].x)):
im = ax.tricontourf(triangulations[1].x, triangulations[1].y, data_2, colors='w')
elif (len(data_2) == len(triangulations[1].xtricenter)):
im = ax.tricontourf(triangulations[1].xtricenter, triangulations[1].ytricenter, data_2, colors='w')
if (len(data_1) == len(triangulations[0].x)):
im = ax.tricontourf(triangulations[0].x, triangulations[0].y, data_1, colors='w')
elif (len(data_1) == len(triangulations[0].xtricenter)):
im = ax.tricontourf(triangulations[0].xtricenter, triangulations[0].ytricenter, data_1, colors='w')
for i in range(len(labels)):
im.collections[i].set_label(labels[i])
im.ax.tick_params(direction='out')
else:
if (len(data_2) == len(triangulations[1].x)):
ax.tricontour(triangulations[1].x, triangulations[1].y, data_2, levels=newlevels, \
linewidths=3, linestyles='dashed', cmap='jet')
elif (len(data_2) == len(triangulations[1].xtricenter)):
ax.tricontour(triangulations[1].xtricenter, triangulations[1].ytricenter, data_2, levels=newlevels, \
linewidths=3, linestyles='dashed', cmap='jet')
if (len(data_1) == len(triangulations[0].x)):
im = ax.tricontour(triangulations[0].x, triangulations[0].y, data_1, levels=newlevels, \
linewidths=3, linestyles='solid', cmap='jet')
elif (len(data_1) == len(triangulations[0].xtricenter)):
im = ax.tricontour(triangulations[0].xtricenter, triangulations[0].ytricenter, data_1, \
levels=newlevels, linewidths=3, linestyles='solid', cmap='jet')
im.set_clim(min(levels),max(levels))
im.ax.tick_params(direction='out')
col = matplotlib.cm.get_cmap('jet')
codecolors = []
newcollection = []
for i in range(len(newlevels)):
thelevelrenorm = (newlevels[i]-min(newlevels))/(max(newlevels)-min(newlevels))
codecolor = col(thelevelrenorm)
codecolors.append(list(codecolor[0:3]))
thepatch = matplotlib.lines.Line2D([], [], color=codecolor, linestyle='solid', linewidth=3, label=labels[i])
newcollection.append(thepatch)
if (len(im.collections)+2 == len(newlevels)):
for i in range(len(im.collections)):
im.collections[i].set_label(labels[i+1])
elif (len(im.collections) == len(newlevels)):
for i in range(len(labels)):
im.collections[i].set_label(labels[i])
#for i in range(len(labels)):
# im.collections[i].set_label(labels[i])
xnodes_1 = triangulations[0].x
ynodes_1 = triangulations[0].y
xmin_1 = min(xnodes_1)
xmax_1 = max(xnodes_1)
ymin_1 = min(ynodes_1)
ymax_1 = max(ynodes_1)
xnodes_2 = triangulations[1].x
ynodes_2 = triangulations[1].y
xmin_2 = min(xnodes_2)
xmax_2 = max(xnodes_2)
ymin_2 = min(ynodes_2)
ymax_2 = max(ynodes_2)
if (geom_kind == 'DIHEDRON'):
triangles_1 = triangulations[0].triangles
themin = ymax_1
for i in range(len(xnodes_1)):
if (abs(xnodes_1[i]-xmin_1)/(xmax_1-xmin_1) < 1e-8):
if (ynodes_1[i] <= themin):
yA1_index = i
themin = ynodes_1[i]
yA1_1 = ynodes_1[yA1_index]
themin = xmax_1
for i in range(len(ynodes_1)):
if (abs(ynodes_1[i]-ymin_1)/(ymax_1-ymin_1) < 1e-8):
if (xnodes_1[i] <= themin):
xA2_index = i
themin = xnodes_1[i]
xA2_1 = xnodes_1[xA2_index]
xslopemin = xA2_1
yslopemin = ymax_1
for i in range(len(triangles_1)):
if (xA2_index in triangles[i]):
x1 = xnodes_1[triangles_1[i][0]]
y1 = ynodes_1[triangles_1[i][0]]
x2 = xnodes_1[triangles_1[i][1]]
y2 = ynodes_1[triangles_1[i][1]]
x3 = xnodes_1[triangles_1[i][2]]
y3 = ynodes_1[triangles_1[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles_1[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles_1[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles_1[i][2]
xA1_1 = (yA1_1-ymin_1)*(xslopemin-xA2_1)/(yslopemin-ymin) + xA2_1
triangles_2 = triangulations[1].triangles
themin = ymax_2
for i in range(len(xnodes_2)):
if (abs(xnodes_2[i]-xmin_2)/(xmax_2-xmin_2) < 1e-8):
if (ynodes_2[i] <= themin):
yA1_index = i
themin = ynodes_2[i]
yA1_2 = ynodes_2[yA1_index]
themin = xmax_2
for i in range(len(ynodes_2)):
if (abs(ynodes_2[i]-ymin_2)/(ymax_2-ymin_2) < 1e-8):
if (xnodes_2[i] <= themin):
xA2_index = i
themin = xnodes_2[i]
xA2_2 = xnodes_2[xA2_index]
xslopemin = xA2_2
yslopemin = ymax_2
for i in range(len(triangles_2)):
if (xA2_index in triangles_2[i]):
x1 = xnodes_2[triangles_2[i][0]]
y1 = ynodes_2[triangles_2[i][0]]
x2 = xnodes_2[triangles_2[i][1]]
y2 = ynodes_2[triangles_2[i][1]]
x3 = xnodes_2[triangles_2[i][2]]
y3 = ynodes_2[triangles_2[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles_2[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles_2[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles_2[i][2]
xA1_2 = (yA1_2-ymin_2)*(xslopemin-xA2_2)/(yslopemin-ymin_2) + xA2_2
if (max([abs(ymin_1-ymin_2), \
abs(xA1_1-xA1_2), abs(yA1_1-yA1_2), abs(xA2_1-xA2_2)]) > 1e-8):
print('Error: domain geometries do not coincide')
exit()
else:
xmin = min([xmin_1, xmin_2])
ymin = ymin_1
xA1 = xA1_1
yA1 = yA1_1
xA2 = xA2_1
xs = [xmin, xA1, xA2, xmin]
ys = [yA1, yA1, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
elif (geom_kind == 'STEP'):
themin = ymax_1
for i in range(len(xnodes_1)):
if (abs(xnodes_1[i]-xmin_1)/(xmax_1-xmin_1) < 1e-8):
if (ynodes_1[i] <= themin):
yminleft_index = i
themin = ynodes_1[i]
yminleft_1 = ynodes_1[yminleft_index]
themin = ymax_1
for i in range(len(xnodes_1)):
if (abs(xnodes_1[i]-xmax_1)/(xmax_1-xmin_1) < 1e-8):
if (ynodes_1[i] <= themin):
yminright_index = i
themin = ynodes_1[i]
yminright_1 = ynodes_1[yminright_index]
themin = ymax_2
for i in range(len(xnodes_2)):
if (abs(xnodes_2[i]-xmin_2)/(xmax_2-xmin_2) < 1e-8):
if (ynodes_2[i] <= themin):
yminleft_index = i
themin = ynodes_2[i]
yminleft_2 = ynodes_2[yminleft_index]
themin = ymax_2
for i in range(len(xnodes_2)):
if (abs(xnodes_2[i]-xmax_2)/(xmax_2-xmin_2) < 1e-8):
if (ynodes_2[i] <= themin):
yminright_index = i
themin = ynodes_2[i]
yminright_2 = ynodes_2[yminright_index]
if (abs(yminleft_1-ymin_1)/(ymax_1-ymin_1) < 1e-8):
ystep_1 = yminright_1
# Step is on the right part of the domain
themax = xmin_1
for i in range(len(ynodes_1)):
if (abs(ynodes_1[i]-ymin_1)/(ymax_1-ymin_1) < 1e-8):
if (xnodes_1[i] >= themax):
xstep_index = i
themax = xnodes_1[i]
xstep_1 = xnodes_1[xstep_index]
if (abs(yminleft_2-ymin_2)/(ymax_2-ymin_2) < 1e-8):
ystep_2 = yminright_2
# Step is on the right part of the domain
themax = xmin_2
for i in range(len(ynodes_2)):
if (abs(ynodes_2[i]-ymin_2)/(ymax_2-ymin_2) < 1e-8):
if (xnodes_2[i] >= themax):
xstep_index = i
themax = xnodes_2[i]
xstep_2 = xnodes_2[xstep_index]
else:
print('Error: Step must be located either of left part or on right part of the domain for both results')
exit()
if (max([abs(xstep_1-xstep_2), abs(ystep_1-ystep_2)]) < 1e-8):
xmax = max([xmax_1, xmax_2])
ymin = min([ymin_1, ymin_2])
xstep = xstep_1
ystep = ystep_1
xs = [xstep, xmax, xmax, xstep]
ys = [ystep, ystep, ymin, ymin]
elif (abs(yminright_1-ymin_1)/(ymax_1-ymin_1) < 1e-8):
ystep_1 = yminleft_1
# Step is on the left part of the domain
themin = xmax_1
for i in range(len(ynodes_1)):
if (abs(ynodes_1[i]-ymin_1)/(ymax_1-ymin_1) < 1e-8):
if (xnodes_1[i] <= themin):
xstep_index = i
themin = xnodes_1[i]
xstep_1 = xnodes_1[xstep_index]
if (abs(yminright_2-ymin_2)/(ymax_2-ymin_2) < 1e-8):
ystep_2 = yminleft_2
# Step is on the left part of the domain
themin = xmax_2
for i in range(len(ynodes_2)):
if (abs(ynodes_2[i]-ymin_2)/(ymax_2-ymin_2) < 1e-8):
if (xnodes_2[i] <= themin):
xstep_index = i
themin = xnodes_2[i]
xstep_2 = xnodes_2[xstep_index]
else:
print('Error: Step must be located either of left part or on right part of the domain for both results')
exit()
if (max([abs(xstep_1-xstep_2), abs(ystep_1-ystep_2)]) < 1e-8):
xmin = min([xmin_1, xmin_2])
ymin = min([ymin_1, ymin_2])
xstep = xstep_1
ystep = ystep_1
xs = [xmin, xstep, xstep, xmin]
ys = [ystep, ystep, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
plt.xlabel("x", fontsize=fs)
plt.ylabel("y", fontsize=fs)
plt.title(aliasfield+" at time t = {:f}".format(time)+" (Re = {:f})".format(Re), fontsize=fs)
ttl = ax.title
ttl.set_position([.5, 1.1])
[oldhandles, labels] = im.ax.get_legend_handles_labels()
if (not(constant)):
oldhandles = newcollection
separator_patch = matplotlib.lines.Line2D([], [], color='w', label='')
leg_1_patch = matplotlib.lines.Line2D([], [], color='black', linestyle='solid', \
linewidth=3, label=aliasfield+' ('+list_alias_diff[0]+')')
leg_2_patch = matplotlib.lines.Line2D([], [], color='black', linestyle='dashed', \
linewidth=3, label=aliasfield+' ('+list_alias_diff[1]+')')
if (min_data_1 < 0):
min_1_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e}'.format(min_data_1))
else:
min_1_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e} '.format(min_data_1))
if (max_data_1 < 0):
max_1_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e}'.format(max_data_1))
else:
max_1_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e} '.format(max_data_1))
if (min_data_2 < 0):
min_2_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e}'.format(min_data_2))
else:
min_2_patch = matplotlib.patches.Patch(color='w', label='Min = {:.5e} '.format(min_data_2))
if (max_data_2 < 0):
max_2_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e}'.format(max_data_2))
else:
max_2_patch = matplotlib.patches.Patch(color='w', label='Max = {:.5e} '.format(max_data_2))
oldhandles.append(separator_patch)
oldhandles.append(leg_1_patch)
oldhandles.append(max_1_patch)
oldhandles.append(min_1_patch)
oldhandles.append(leg_2_patch)
oldhandles.append(max_2_patch)
oldhandles.append(min_2_patch)
lgd = plt.legend(handles=oldhandles, fontsize=fs, ncol=1, loc='upper left', \
frameon=False, bbox_to_anchor=tuple_anchor_legend)
if (len(list_bbox) == 4):
if (list_bbox[0] == 'default'):
xminset = min([xmin_1, xmin_2])
else:
xminset = list_bbox[0]
if (list_bbox[1] == 'default'):
xmaxset = max([xmax_1, xmax_2])
else:
xmaxset = list_bbox[1]
if (list_bbox[2] == 'default'):
yminset = min([ymin_1, ymin_2])
else:
yminset = list_bbox[2]
if (list_bbox[3] == 'default'):
ymaxset = max([ymax_1, ymax_2])
else:
ymaxset = list_bbox[3]
ax.set_xlim([xminset, xmaxset])
ax.set_ylim([yminset, ymaxset])
else:
print('Error: list_bbox must be a list of size 4')
exit()
plt.draw()
if (not(outputpng == '')):
fig.savefig(outputpng, bbox_extra_artists=(lgd,ttl,), bbox_inches='tight')
del data_1, data_2, newlevels, xnodes_1, ynodes_1, xnodes_2, ynodes_2
if (geom_kind == 'DIHEDRON'):
del triangles_1, triangles_2#, x, y1, y2
fig.clf()
plt.close()
h5file_1.close()
h5file_2.close()
def plot_two_contour(inputh5, field, list_alias_diff, \
aliasfield='', \
levels=[], \
list_bbox=['default']*4, \
tuple_anchor_legend=(1.01,1.01), \
tuple_inches_figsize=(12.,12.), \
outputpng='./temp.png', \
forcesave=False):
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
f1 = h5py.File(inputh5[0], 'r')
f2 = h5py.File(inputh5[1], 'r')
meshh51 = f1.attrs['Mesh'].decode('UTF-8')
meshh52 = f2.attrs['Mesh'].decode('UTF-8')
if (field in list_Pp):
method1 = f1.attrs['Method_pressure'].decode('UTF-8')
method2 = f2.attrs['Method_pressure'].decode('UTF-8')
elif (field in list_Pv):
method1 = f1.attrs['Method_velocity'].decode('UTF-8')
method2 = f2.attrs['Method_velocity'].decode('UTF-8')
elif (field in list_Prho):
method1 = f1.attrs['Method_density'].decode('UTF-8')
method2 = f2.attrs['Method_density'].decode('UTF-8')
f1.close()
f2.close()
if (len(levels) == 0):
# Need to provide isovalues
theminmax_1 = extract_minmax(inputh5[0], field)
theminmax_2 = extract_minmax(inputh5[1], field)
newlevels = numpy.linspace(min([min(theminmax_1), min(theminmax_2)]), \
max([max(theminmax_1), max(theminmax_2)]), 10, endpoint=True)
else:
newlevels = levels.copy()
if (aliasfield == ''):
newaliasfield = field
else:
newaliasfield = aliasfield
# We verify if the mesh file is correct
if (not(os.path.isfile(meshh51))):
# The path to the mesh file may be wrong
testmeshh51 = inputh5[0].split(inputh5[0].split('/')[-1])[0] + meshh51.split('/')[-1]
if (not(os.path.isfile(testmeshh51))):
print('Problem: cannot find the mesh file associated to '+inputh5[0])
exit()
else:
# Mesh file found
meshh51 = testmeshh51
if (not(os.path.isfile(meshh52))):
# The path to the mesh file may be wrong
testmeshh52 = inputh5[1].split(inputh5[1].split('/')[-1])[0] + meshh52.split('/')[-1]
if (not(os.path.isfile(testmeshh52))):
print('Problem: cannot find the mesh file associated to '+inputh5[1])
exit()
else:
# Mesh file found
meshh52 = testmeshh52
# Load mesh
triangulation1 = build_triangulation(meshh51, method1)
triangulation2 = build_triangulation(meshh52, method2)
triangulations = [triangulation1, triangulation2]
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with file " + outputpng + ": the file already exists. Please choose another name")
exit()
plot_two_contour_withtri(inputh5, field, newaliasfield, list_alias_diff, newlevels, list_bbox, \
triangulations, tuple_anchor_legend, tuple_inches_figsize, outputpng)
def plot_contour_set(inputh5set, field, \
aliasfield='', \
nlevels=10, \
fix_isovalues=True, \
list_bbox=['default']*4, \
tuple_anchor_legend=(1.01,1.01), \
tuple_inches_figsize=(12., 12.), \
outputdir='.', \
forcesave=False):
ndigits = int(numpy.ceil(numpy.log10(len(inputh5set)))+1)
if (not(os.path.isdir(outputdir))):
os.mkdir(outputdir)
print("Directory "+outputdir+" created")
creation = True
else:
print("Directory "+outputdir+" already exists")
creation = False
print(str(len(inputh5set)) + ' PNG files have to be generated.')
if (aliasfield == ''):
newaliasfield = field
else:
newaliasfield = aliasfield
if (fix_isovalues):
# Compute reference isovalues that will be the same for the whole file set
# Need to find the min and max values of the field over the whole file set
list_min = []
list_max = []
for inputh5 in inputh5set:
theminmax = extract_minmax(inputh5, field)
list_min.append(theminmax[0])
list_max.append(theminmax[1])
list_minmax = [min(list_min), max(list_max)]
# Use these values to build the isovalues
levelsref = numpy.linspace(min(list_minmax), max(list_minmax), nlevels, endpoint=True)
else:
# The isovalues depend on the opened file
levelsref = []
# if (len(list_minmax) == 2):
# if ((max(list_minmax)-min(list_minmax)) > 1e-12*max([abs(max(list_minmax)),abs(min(list_minmax))])):
# levelsref = numpy.linspace(min(list_minmax), max(list_minmax), nlevels, endpoint=True)
# else:
# levelsref = []
# else:
# levelsref = []
meshh5 = ''
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
for inputh5 in inputh5set:
if (levelsref == []):
# No reference isovalues: we build them from the current file
[minval, maxval] = extract_minmax(inputh5, field)
if (maxval-minval > 1e-12*max([abs(maxval),abs(minval)])):
levels = numpy.linspace(minval, maxval, nlevels, endpoint=True)
else:
levels = []
else:
# Use reference isovalues
levels = levelsref.copy()
# if required, load a triangulation
f = h5py.File(inputh5, 'r')
newmeshh5 = f.attrs['Mesh'].decode('UTF-8')
if (field in list_Pp):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif (field in list_Pv):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif (field in list_Prho):
method = f.attrs['Method_density'].decode('UTF-8')
f.close()
# We verify if the mesh file is correct
if (not(os.path.isfile(newmeshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + newmeshh5.split('/')[-1]
#print(testmeshh5)
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
newmeshh5 = testmeshh5
if (not(newmeshh5 == meshh5)):
print('load new triangulation')
# Need to load a new mesh
triangulation = build_triangulation(newmeshh5, method)
meshh5 = newmeshh5
# At this point, a triangulation should have been loaded
# Find the prefix of the namefile
trunch5 = inputh5.split('/')[-1]
extensionh5 = trunch5.split('_')[-1]
prefixh5 = trunch5.split(extensionh5)[0][:-1]
indexh5 = int(extensionh5.split('.')[0])
outputpng = outputdir + '/' + prefixh5 + '_' + str(indexh5).rjust(ndigits, '0') + '_' + \
aliasfield + '.png'
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with directory " + outputdir + ": the file " + \
outputpng + " already exists (and probably other file). Please choose another directory")
exit()
plot_contour_withtri(inputh5, field, newaliasfield, levels, list_bbox, triangulation, \
tuple_anchor_legend, tuple_inches_figsize, outputpng)
print('File '+outputpng+' saved')
del levels
del levelsref, triangulation
def plot_pseudocolor_withtri(inputh5, field, aliasfield, list_minmax, list_bbox, triangulation, \
tuple_leg_position, tuple_inches_figsize, outputpng):
h5file = h5py.File(inputh5, 'r')
geom_kind = h5file.attrs['Geometry'].decode('UTF-8')
if (field in h5file.keys()):
data = h5file[field].value
else:
print('Error: cannot find the field '+field+' in file '+inputh5)
exit()
if ((not(data.shape[0] == 1)) and (data.shape[1] == 1)):
# Need to transpose the data
data = numpy.reshape(data, (1,data.shape[0]))[0]
else:
data = data[0]
min_data = min(data)
max_data = max(data)
if (len(list_minmax) == 2):
min_data_set = min(list_minmax)
max_data_set = max(list_minmax)
else:
min_data_set = min_data
max_data_set = max_data
if ('Time' in h5file.keys()):
time = h5file['Time'].value[0]
else:
print('Error: cannot find the field ''Time'' in file '+inputh5)
exit()
if ('Reynolds_number' in h5file.keys()):
Re = h5file['Reynolds_number'].value[0]
else:
print('Error: cannot find the field ''Reynolds_number'' in file '+inputh5)
exit()
if ((max_data-min_data) > 1e-12*max([abs(max_data),abs(min_data)])):
constant = False
levels = numpy.linspace(min_data_set, max_data_set, 256, endpoint=True)
newlevels = levels.copy()
newlevels[0] += (levels[-1]-levels[0])*0.000001
newlevels[-1] -= (levels[-1]-levels[0])*0.000001
labels = []
else:
constant = True
labels = [aliasfield+" = {:.5e}".format(min_data), "(Constant)"]
fs = 20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
matplotlib.rcParams['font.family'] = 'monospace'
fig = plt.figure(figsize=tuple_inches_figsize)
ax = plt.gca()
ax.set_aspect('equal')
if (constant):
if (len(data) == len(triangulation.x)):
im = ax.tricontourf(triangulation.x, triangulation.y, data, colors='greenyellow')
elif (len(data) == len(triangulation.xtricenter)):
im = ax.tricontourf(triangulation.xtricenter, triangulation.ytricenter, data, colors='greenyellow')
im.ax.tick_params(direction='out')
cb = fig.colorbar(im, aspect=1, format=aliasfield+" = %.5e\n(Constant)", \
ticks=[min_data])
cb.ax.tick_params(axis='y', direction='in', color='greenyellow')
cb.ax.set_position(pos=(tuple_leg_position[0], tuple_leg_position[1], 0.05, 0.05))
else:
if (len(data) == len(triangulation.x)):
im = ax.tricontourf(triangulation.x, triangulation.y, data, levels=newlevels, cmap='jet')
elif (len(data) == len(triangulation.xtricenter)):
im = ax.tricontourf(triangulation.xtricenter, triangulation.ytricenter, data, \
levels=newlevels, cmap='jet')
im.ax.tick_params(direction='out')
cb = fig.colorbar(im, aspect=10, format="%.5e", \
ticks=numpy.linspace(newlevels[0],newlevels[-1],7,endpoint=True).tolist())
cb.ax.tick_params(axis='y', direction='out')
cb.ax.set_position(pos=tuple_leg_position)
xnodes = triangulation.x
ynodes = triangulation.y
xmin = min(xnodes)
xmax = max(xnodes)
ymin = min(ynodes)
ymax = max(ynodes)
if (geom_kind == 'DIHEDRON'):
triangles = triangulation.triangles
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yA1_index = i
themin = ynodes[i]
yA1 = ynodes[yA1_index]
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xA2_index = i
themin = xnodes[i]
xA2 = xnodes[xA2_index]
xslopemin = xA2
yslopemin = ymax
for i in range(len(triangles)):
if (xA2_index in triangles[i]):
x1 = xnodes[triangles[i][0]]
y1 = ynodes[triangles[i][0]]
x2 = xnodes[triangles[i][1]]
y2 = ynodes[triangles[i][1]]
x3 = xnodes[triangles[i][2]]
y3 = ynodes[triangles[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles[i][2]
xA1 = (yA1-ymin)*(xslopemin-xA2)/(yslopemin-ymin) + xA2
xs = [xmin, xA1, xA2, xmin]
ys = [yA1, yA1, ymin, ymin]
ax.fill(xs, ys, 'white')
elif (geom_kind == 'STEP'):
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminleft_index = i
themin = ynodes[i]
yminleft = ynodes[yminleft_index]
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmax)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminright_index = i
themin = ynodes[i]
yminright = ynodes[yminright_index]
if (abs(yminleft-ymin)/(ymax-ymin) < 1e-8):
ystep = yminright
# Step is on the right part of the domain
themax = xmin
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] >= themax):
xstep_index = i
themax = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xstep, xmax, xmax, xstep]
ys = [ystep, ystep, ymin, ymin]
elif (abs(yminright-ymin)/(ymax-ymin) < 1e-8):
ystep = yminleft
# Step is on the left part of the domain
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xstep_index = i
themin = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xmin, xstep, xstep, xmin]
ys = [ystep, ystep, ymin, ymin]
ax.fill(xs, ys, 'white')
plt.xlabel("x", fontsize=fs)
plt.ylabel("y", fontsize=fs)
plt.title(aliasfield+" at time t = {:f}".format(time)+" (Re = {:f})".format(Re), fontsize=fs)
ttl = ax.title
ttl.set_position([.5, 1.1])
tminmax = ''
tminmax = tminmax+'Min('+aliasfield+') = {:.5e}'.format(min_data)
if (min_data >= 0):
tminmax = tminmax+' '
tminmax = tminmax+'\n'
tminmax = tminmax+'Max('+aliasfield+') = {:.5e}'.format(max_data)
if (max_data >= 0):
tminmax = tminmax+' '
[oldhandles, labels] = cb.ax.get_legend_handles_labels()
lgd = cb.ax.legend(handles=oldhandles, fontsize=fs, ncol=1, loc='upper left', frameon=False, \
bbox_to_anchor=(0,0))
if (constant):
tt = cb.ax.text(0, -2, tminmax, fontsize=fs, linespacing=1.5)
else:
tt = cb.ax.text(0, -0.25, tminmax, fontsize=fs, linespacing=1.5)
if (len(list_bbox) == 4):
if (list_bbox[0] == 'default'):
xminset = xmin
else:
xminset = list_bbox[0]
if (list_bbox[1] == 'default'):
xmaxset = xmax
else:
xmaxset = list_bbox[1]
if (list_bbox[2] == 'default'):
yminset = ymin
else:
yminset = list_bbox[2]
if (list_bbox[3] == 'default'):
ymaxset = ymax
else:
ymaxset = list_bbox[3]
ax.set_xlim([xminset, xmaxset])
ax.set_ylim([yminset, ymaxset])
else:
print('Error: list_bbox must be a list of size 4')
exit()
plt.draw()
if (not(outputpng == '')):
fig.savefig(outputpng, bbox_extra_artists=(lgd,ttl,tt,), bbox_inches='tight')
del data
if (not(constant)):
del newlevels
fig.clf()
plt.close()
h5file.close()
def plot_pseudocolor(inputh5, field, aliasfield='', list_minmax=[], \
list_bbox=['default']*4, \
tuple_leg_position=(0.8,0.45,0.2,0.3), \
tuple_inches_figsize=(12.,12.), \
outputpng='./temp.png', \
forcesave=False):
if (aliasfield == ''):
newaliasfield = field
else:
newaliasfield = aliasfield
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
f = h5py.File(inputh5, 'r')
meshh5 = f.attrs['Mesh'].decode('UTF-8')
if (field in list_Pp):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif (field in list_Pv):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif (field in list_Prho):
method = f.attrs['Method_density'].decode('UTF-8')
f.close()
# We verify if the mesh file is correct
if (not(os.path.isfile(meshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + meshh5.split('/')[-1]
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
meshh5 = testmeshh5
# Load mesh
triangulation = build_triangulation(meshh5, method)
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with file " + outputpng + ": the file already exists. Please choose another name")
exit()
plot_pseudocolor_withtri(inputh5, field, newaliasfield, list_minmax, list_bbox, triangulation, \
tuple_leg_position, tuple_inches_figsize, outputpng)
def plot_pseudocolor_set(inputh5set, field, \
aliasfield='', \
list_minmax=[], \
list_bbox=['default']*4, \
tuple_leg_position=(0.8,0.45,0.2,0.3), \
tuple_inches_figsize=(12., 12.), \
outputdir='.', \
forcesave=False):
if (aliasfield == ''):
newaliasfield = field
else:
newaliasfield = aliasfield
print(str(len(inputh5set)) + ' PNG files have to be generated.')
ndigits = int(numpy.ceil(numpy.log10(len(inputh5set)))+1)
if (not(os.path.isdir(outputdir))):
os.mkdir(outputdir)
print("Directory "+outputdir+" created")
creation = True
else:
print("Directory "+outputdir+" already exists")
creation = False
if (len(list_minmax) == 0):
# Compute min and max bound that will be the same for the whole file set
list_min = []
list_max = []
for inputh5 in inputh5set:
theminmax = extract_minmax(inputh5, field)
list_min.append(theminmax[0])
list_max.append(theminmax[1])
newlist_minmax = [min(list_min), max(list_max)]
elif (len(list_minmax) == 2):
newlist_minmax = [min(list_minmax), max(list_minmax)]
else:
print('Problem with the argument list_minmax in plot_pseudocolor_set: it should be of size 0 or 2')
exit()
meshh5 = ''
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
for inputh5 in inputh5set:
# if required, load a triangulation
f = h5py.File(inputh5, 'r')
newmeshh5 = f.attrs['Mesh'].decode('UTF-8')
if (field in list_Pp):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif (field in list_Pv):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif (field in list_Prho):
method = f.attrs['Method_density'].decode('UTF-8')
f.close()
# We verify if the mesh file is correct
if (not(os.path.isfile(newmeshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + newmeshh5.split('/')[-1]
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
newmeshh5 = testmeshh5
if (not(newmeshh5 == meshh5)):
print('load new triangulation')
# Need to load a new mesh
triangulation = build_triangulation(newmeshh5, method)
meshh5 = newmeshh5
# At this point, a triangulation should have been loaded
# Find the prefix of the namefile
trunch5 = inputh5.split('/')[-1]
extensionh5 = trunch5.split('_')[-1]
prefixh5 = trunch5.split(extensionh5)[0][:-1]
indexh5 = int(extensionh5.split('.')[0])
outputpng = outputdir + '/' + prefixh5 + '_' + str(indexh5).rjust(ndigits, '0') + '_' + \
aliasfield + '.png'
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with directory " + outputdir + ": the file " + \
outputpng + " already exists (and probably other file). Please choose another directory")
exit()
plot_pseudocolor_withtri(inputh5, field, newaliasfield, newlist_minmax, list_bbox, triangulation, \
tuple_leg_position, tuple_inches_figsize, outputpng)
print('File '+outputpng+' saved')
del triangulation
def plot_vectorfield_withtri(inputh5, field_x, field_y, narrows, scale, width, list_bbox, \
triangulation, tuple_leg_position, tuple_inches_figsize, outputpng, aliasfield=''):
if (aliasfield == ''):
aliasx = field_x
aliasy = field_y
alias = '('+field_x+','+field_y+')'
else:
aliasx = aliasfield+'_x'
aliasy = aliasfield+'_y'
alias = aliasfield
h5file = h5py.File(inputh5, 'r')
geom_kind = h5file.attrs['Geometry'].decode('UTF-8')
if (field_x in h5file.keys()):
data_x = h5file[field_x].value
else:
print('Error: cannot find the field '+field_x+' in file '+inputh5)
exit()
data_x = numpy.asarray(data_x)
if ((not(data_x.shape[0] == 1)) and (data_x.shape[1] == 1)):
# Need to transpose the data
data_x = numpy.reshape(data_x, (1,data_x.shape[0]))[0]
else:
data_x = data_x[0]
if (field_y in h5file.keys()):
data_y = h5file[field_y].value
else:
print('Error: cannot find the field '+field_y+' in file '+inputh5)
exit()
data_y = h5file[field_y].value
if ((not(data_y.shape[0] == 1)) and (data_y.shape[1] == 1)):
# Need to transpose the data
data_y = numpy.reshape(data_y, (1,data_y.shape[0]))[0]
else:
data_y = data_y[0]
if (len(data_x) == len(data_y)):
if (len(data_x) == len(triangulation.x)):
Xnodes = triangulation.x
Ynodes = triangulation.y
elif (len(data_x) == len(triangulation.xtricenter)):
Xnodes = triangulation.xtricenter
Ynodes = triangulation.ytricenter
else:
print('Error: cannot conclude if the data is cell-centered or node-centered')
exit()
else:
print('Error: fields '+field_x+' and '+field_y+' must have same size')
min_data_x = min(data_x)
max_data_x = max(data_x)
min_data_y = min(data_y)
max_data_y = max(data_y)
xmin = min(triangulation.x)
xmax = max(triangulation.x)
ymin = min(triangulation.y)
ymax = max(triangulation.y)
bb = numpy.asarray([[xmin,ymin], [xmax,ymax]])
constant = not((abs(max_data_x-min_data_x) > 1e-12*max([abs(max_data_x),abs(min_data_x)])) or \
(abs(max_data_y-min_data_y) > 1e-12*max([abs(max_data_y),abs(min_data_y)])))
data_color = []
for i in range(len(data_x)):
data_color.append(numpy.sqrt(data_x[i]*data_x[i]+data_y[i]*data_y[i]))
if ('Time' in h5file.keys()):
time = h5file['Time'].value[0]
else:
print('Error: cannot find the field ''Time'' in file '+inputh5)
if ('Reynolds_number' in h5file.keys()):
Re = h5file['Reynolds_number'].value[0]
else:
print('Error: cannot find the field ''Reynolds_number'' in file '+inputh5)
fs = 20
plt.rc('xtick',labelsize=fs)
plt.rc('ytick',labelsize=fs)
matplotlib.rcParams['font.family'] = 'monospace'
fig = plt.figure(figsize=tuple_inches_figsize)
ax = plt.gca()
ax.set_aspect('equal')
step = len(data_x)//narrows
Cmin_ind = numpy.argmin(data_color).item()
Cmax_ind = numpy.argmax(data_color).item()
tab = Xnodes.tolist()[0:len(Xnodes):step]
tab.append(Xnodes[Cmin_ind])
tab.append(Xnodes[Cmax_ind])
X = numpy.asarray(tab)
tab = Ynodes.tolist()[0:len(Ynodes):step]
tab.append(Ynodes[Cmin_ind])
tab.append(Ynodes[Cmax_ind])
Y = numpy.asarray(tab)
tab = data_x.tolist()[0:len(data_x):step]
tab.append(data_x[Cmin_ind])
tab.append(data_x[Cmax_ind])
VX = numpy.asarray(tab)
tab = data_y.tolist()[0:len(data_y):step]
tab.append(data_y[Cmin_ind])
tab.append(data_y[Cmax_ind])
VY = numpy.asarray(tab)
tab = data_color[0:len(data_color):step]
tab.append(data_color[Cmin_ind])
tab.append(data_color[Cmax_ind])
C = numpy.asarray(tab)
if (not(constant)):
im = ax.quiver(X, Y, VX, VY, C, units='xy', scale=scale, zorder=3, \
cmap='jet', width=width, headwidth=10, headlength=10., headaxislength=6)
im.ax.set_xlim(xmin,xmax)
im.ax.set_ylim(ymin,ymax)
im.ax.tick_params(direction='out')
cb = fig.colorbar(im, aspect=10, format="%.5e", \
ticks=numpy.linspace(min(data_color),max(data_color),7,endpoint=True).tolist())
cb.ax.tick_params(axis='y', direction='out')
cb.ax.set_position(pos=tuple_leg_position)
else:
#im = ax.quiver(X, Y, VX, VY, C, units='xy', scale=scale, zorder=3, \
# cmap='jet', width=width, headwidth=10, headlength=10., headaxislength=6)
im = ax.tricontourf(triangulation.x, triangulation.y, data_color, colors='w')
im.ax.tick_params(direction='out')
cb = fig.colorbar(im, aspect=1, format="|"+alias+"| = %.5e\n(Constant)", \
ticks=[data_color[0]])
cb.ax.tick_params(axis='y', direction='in', color='w')
cb.ax.set_position(pos=(tuple_leg_position[0], tuple_leg_position[1], 0.05, 0.05))
xnodes = triangulation.x
ynodes = triangulation.y
xmin = min(xnodes)
xmax = max(xnodes)
ymin = min(ynodes)
ymax = max(ynodes)
if (geom_kind == 'DIHEDRON'):
triangles = triangulation.triangles
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yA1_index = i
themin = ynodes[i]
yA1 = ynodes[yA1_index]
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xA2_index = i
themin = xnodes[i]
xA2 = xnodes[xA2_index]
xslopemin = xA2
yslopemin = ymax
for i in range(len(triangles)):
if (xA2_index in triangles[i]):
x1 = xnodes[triangles[i][0]]
y1 = ynodes[triangles[i][0]]
x2 = xnodes[triangles[i][1]]
y2 = ynodes[triangles[i][1]]
x3 = xnodes[triangles[i][2]]
y3 = ynodes[triangles[i][2]]
if ((x1 < xslopemin) and (y1 < yslopemin)):
xslopemin = x1
yslopemin = y1
slopemin_index = triangles[i][0]
if ((x2 < xslopemin) and (y2 < yslopemin)):
xslopemin = x2
yslopemin = y2
slopemin_index = triangles[i][1]
if ((x3 < xslopemin) and (y3 < yslopemin)):
xslopemin = x3
yslopemin = y3
slopemin_index = triangles[i][2]
xA1 = (yA1-ymin)*(xslopemin-xA2)/(yslopemin-ymin) + xA2
xs = [xmin, xA1, xA2, xmin]
ys = [yA1, yA1, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
elif (geom_kind == 'STEP'):
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmin)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminleft_index = i
themin = ynodes[i]
yminleft = ynodes[yminleft_index]
themin = ymax
for i in range(len(xnodes)):
if (abs(xnodes[i]-xmax)/(xmax-xmin) < 1e-8):
if (ynodes[i] <= themin):
yminright_index = i
themin = ynodes[i]
yminright = ynodes[yminright_index]
if (abs(yminleft-ymin)/(ymax-ymin) < 1e-8):
ystep = yminright
# Step is on the right part of the domain
themax = xmin
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] >= themax):
xstep_index = i
themax = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xstep, xmax, xmax, xstep]
ys = [ystep, ystep, ymin, ymin]
elif (abs(yminright-ymin)/(ymax-ymin) < 1e-8):
ystep = yminleft
# Step is on the left part of the domain
themin = xmax
for i in range(len(ynodes)):
if (abs(ynodes[i]-ymin)/(ymax-ymin) < 1e-8):
if (xnodes[i] <= themin):
xstep_index = i
themin = xnodes[i]
xstep = xnodes[xstep_index]
xs = [xmin, xstep, xstep, xmin]
ys = [ystep, ystep, ymin, ymin]
ax.fill(xs, ys, 'lightgray')
plt.xlabel("x", fontsize=fs)
plt.ylabel("y", fontsize=fs)
plt.title(alias+" at time t = {:f}".format(time)+" (Re = {:f})".format(Re), fontsize=fs)
ttl = ax.title
ttl.set_position([.5, 1.1])
tminmax = ''
tminmax = tminmax+'Min('+aliasx+') = {:.5e}'.format(min_data_x)
if (min_data_x >= 0):
tminmax = tminmax+' '
tminmax = tminmax+'\n'
tminmax = tminmax+'Max('+aliasx+') = {:.5e}'.format(max_data_x)
if (max_data_x >= 0):
tminmax = tminmax+' '
tminmax = tminmax+'\n'
tminmax = tminmax+'Min('+aliasy+') = {:.5e}'.format(min_data_y)
if (min_data_y >= 0):
tminmax = tminmax+' '
tminmax = tminmax+'\n'
tminmax = tminmax+'Max('+aliasy+') = {:.5e}'.format(max_data_y)
if (max_data_y >= 0):
tminmax = tminmax+' '
[oldhandles, labels] = cb.ax.get_legend_handles_labels()
lgd = cb.ax.legend(handles=oldhandles, fontsize=fs, ncol=1, loc='upper left', frameon=False, \
bbox_to_anchor=(0,0))
if (constant):
tt = cb.ax.text(0, -3, tminmax, fontsize=fs, linespacing=1.5)
else:
tt = cb.ax.text(0, -0.5, tminmax, fontsize=fs, linespacing=1.5)
if (len(list_bbox) == 4):
if (list_bbox[0] == 'default'):
xminset = xmin
else:
xminset = list_bbox[0]
if (list_bbox[1] == 'default'):
xmaxset = xmax
else:
xmaxset = list_bbox[1]
if (list_bbox[2] == 'default'):
yminset = ymin
else:
yminset = list_bbox[2]
if (list_bbox[3] == 'default'):
ymaxset = ymax
else:
ymaxset = list_bbox[3]
ax.set_xlim([xminset, xmaxset])
ax.set_ylim([yminset, ymaxset])
else:
print('Error: list_bbox must be a list of size 4')
exit()
plt.draw()
if (not(outputpng == '')):
fig.savefig(outputpng, bbox_extra_artists=(lgd,ttl,tt,), bbox_inches='tight')
del data_x, data_y, xnodes, ynodes
if (geom_kind == 'DIHEDRON'):
del triangles
fig.clf()
plt.close()
h5file.close()
def plot_vectorfield(inputh5, field_x, field_y, \
aliasfield='', \
narrows=200, \
scale=2., \
width=0.01, \
list_bbox=['default']*4, \
tuple_leg_position=(0.8,0.45,0.2,0.3), \
tuple_inches_figsize=(12.,12.), \
outputpng='./temp.png', \
forcesave=False):
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
f = h5py.File(inputh5, 'r')
meshh5 = f.attrs['Mesh'].decode('UTF-8')
if ((field_x in list_Pp) and (field_y in list_Pp)):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif ((field_x in list_Pv) and (field_y in list_Pv)):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif ((field_x in list_Prho) and (field_y in list_Prho)):
method = f.attrs['Method_density'].decode('UTF-8')
else:
print('Error: scalar fields '+field_x+' and '+field_y+' must be linked to the same finite element nodes')
exit()
f.close()
# We verify if the mesh file is correct
if (not(os.path.isfile(meshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + meshh5.split('/')[-1]
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
meshh5 = testmeshh5
# Load mesh
triangulation = build_triangulation(meshh5, method)
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with file " + outputpng + ": the file already exists. Please choose another name")
exit()
plot_vectorfield_withtri(inputh5, field_x, field_y, narrows, scale, width, list_bbox, triangulation, \
tuple_leg_position, tuple_inches_figsize, outputpng, aliasfield)
def plot_vectorfield_set(inputh5set, field_x, field_y, \
aliasfield='', \
narrows=200, \
scale=2., \
width=0.01, \
list_bbox=['default']*4, \
tuple_leg_position=(0.8,0.45,0.2,0.3), \
tuple_inches_figsize=(12.,12.), \
outputdir='.', \
forcesave=False):
ndigits = int(numpy.ceil(numpy.log10(len(inputh5set)))+1)
if (not(os.path.isdir(outputdir))):
os.mkdir(outputdir)
print("Directory "+outputdir+" created")
creation = True
else:
print("Directory "+outputdir+" already exists")
creation = False
print(str(len(inputh5set)) + ' PNG files have to be generated.')
meshh5 = ''
list_Pp = ['Pressure', 'Pressure_gradient_x', 'Pressure_gradient_y', 'Streamlines', 'Vorticity', \
'Velocity_x_dx', 'Velocity_x_dy', 'Velocity_y_dx', 'Velocity_y_dy', 'Shear_rate', \
'r', \
'Pressure_exact', 'Pressure_gradient_x_exact', 'Pressure_gradient_y_exact', \
'Streamlines_exact', 'Vorticity_exact', \
'Velocity_x_dx_exact', 'Velocity_x_dy_exact', 'Velocity_y_dx_exact', 'Velocity_y_dy_exact', \
'Shear_rate_exact']
list_Pv = ['Velocity_x', 'Velocity_y', 'Velocity_magnitude', 'wx', 'wy', \
'Velocity_x_exact', 'Velocity_y_exact', 'Velocity_magnitude_exact']
list_Prho = ['Density', 'Density_exact']
for inputh5 in inputh5set:
# if required, load a triangulation
f = h5py.File(inputh5, 'r')
newmeshh5 = f.attrs['Mesh'].decode('UTF-8')
if ((field_x in list_Pp) and (field_y in list_Pp)):
method = f.attrs['Method_pressure'].decode('UTF-8')
elif ((field_x in list_Pv) and (field_y in list_Pv)):
method = f.attrs['Method_velocity'].decode('UTF-8')
elif ((field_x in list_Prho) and (field_y in list_Prho)):
method = f.attrs['Method_density'].decode('UTF-8')
else:
print('Error: scalar fields '+field_x+' and '+field_y+' must be linked to the same finite element nodes')
exit()
f.close()
# We verify if the mesh file is correct
if (not(os.path.isfile(newmeshh5))):
# The path to the mesh file may be wrong
testmeshh5 = inputh5.split(inputh5.split('/')[-1])[0] + newmeshh5.split('/')[-1]
if (not(os.path.isfile(testmeshh5))):
print('Problem: cannot find the mesh file associated to '+inputh5)
exit()
else:
# Mesh file found
newmeshh5 = testmeshh5
if (not(newmeshh5 == meshh5)):
print('load new triangulation')
# Need to load a new mesh
triangulation = build_triangulation(newmeshh5, method)
meshh5 = newmeshh5
# At this point, a triangulation should have been loaded
# Find the prefix of the namefile
trunch5 = inputh5.split('/')[-1]
extensionh5 = trunch5.split('_')[-1]
prefixh5 = trunch5.split(extensionh5)[0][:-1]
indexh5 = int(extensionh5.split('.')[0])
outputpng = outputdir + '/' + prefixh5 + '_' + str(indexh5).rjust(ndigits, '0') + '_' + \
aliasfield + '.png'
if (os.path.isfile(outputpng) and not(forcesave)):
print("Problem with directory " + outputdir + ": the file " + \
outputpng + " already exists (and probably other file). Please choose another directory")
exit()
plot_vectorfield_withtri(inputh5, field_x, field_y, narrows, scale, width, list_bbox, triangulation, \
tuple_leg_position, tuple_inches_figsize, outputpng, aliasfield)
print('File '+outputpng+' saved')
del triangulation
