import os
from issm.expread import expread
import numpy as np
from issm.project2d import project2d
#from issm.InterpFromMesh2d import InterpFromMesh2d
from issm.InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
from issm.InterpFromMeshToMesh3d import InterpFromMeshToMesh3d
[docs]def SectionValues(md,data,infile,resolution):
'''
compute the value of a field on a section
This routine gets the value of a given field of the model on points
given in the file infile (Argus type file). Resolution must be a list
[horizontal_resolution, vertical_resolution]
Usage:
[elements,x,y,z,s,data]=SectionValues(md,data,filename,resolution)
[elements,x,y,z,s,data]=SectionValues(md,data,profile_structure,resolution)
'''
if os.path.isfile(infile):
profile=expread(infile)[0]
nods=profile['nods']
x=profile['x']
y=profile['y']
else:
raise IOError('file %s not found' % infile)
#get the specified resolution
if len(resolution)!=2:
raise ValueError('SectionValues error message: Resolution must be a list [horizontal_resolution, vertical_resolution]')
else:
res_h=resolution[0]
if md.mesh.domaintype().lower() == '3d':
if isinstance(resolution[1],int) or isinstance(resolution[1],float):
res_v=resolution[1]
else:
raise ValueError('SectionValues error: resolution must be a length-2 list of integers or floats')
#initialization
X=np.array([]) #X-coordinate
Y=np.array([]) #Y-coordinate
S=np.array([0.]) #curvilinear coordinate
for i in xrange(nods-1):
x_start=x[i]
x_end=x[i+1]
y_start=y[i]
y_end=y[i+1]
s_start=S[-1]
length_segment=np.sqrt((x_end-x_start)**2+(y_end-y_start)**2)
portion=np.ceil(length_segment/res_h)
x_segment=np.zeros(portion)
y_segment=np.zeros(portion)
s_segment=np.zeros(portion)
for j in xrange(int(portion)):
x_segment[j]=x_start+(j)*(x_end-x_start)/portion
y_segment[j]=y_start+(j)*(y_end-y_start)/portion
s_segment[j]=s_start+j*length_segment/portion
#plug into X and Y
X=np.append(X,x_segment)
Y=np.append(Y,y_segment)
S=np.append(S,s_segment)
X=np.append(X,x[nods-1])
Y=np.append(Y,y[nods-1])
#Number of nodes:
numberofnodes=X.shape[0]
#Compute Z
Z=np.zeros(numberofnodes)
#New mesh and Data interpolation
if '2d' in md.mesh.domaintype().lower():
#Interpolation of data on specified points
#data_interp=InterpFromMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y)[0]
data_interp=InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,data,X,Y)[0]
#data_interp=griddata(md.mesh.x,md.mesh.y,data,X,Y)
#Compute index
index=np.array([range(1,numberofnodes),range(2,numberofnodes+1)]).T
else:
#vertically extrude mesh
#Get base and surface for each 2d point, offset to make sure that it is inside the glacier system
offset=1.e-3
base=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.base,1),X,Y)[0]+offset
base=base.reshape(-1,)
surface=InterpFromMeshToMesh2d(md.mesh.elements2d,md.mesh.x2d,md.mesh.y2d,project2d(md,md.geometry.surface,1),X,Y)[0]-offset
surface=surface.reshape(-1,)
#Some useful parameters
layers=int(np.ceil(np.mean(md.geometry.thickness)/res_v))
nodesperlayer=int(numberofnodes)
nodestot=int(nodesperlayer*layers)
elementsperlayer=int(nodesperlayer-1)
elementstot=int((nodesperlayer-1)*(layers-1))
#initialization
X3=np.zeros(nodesperlayer*layers)
Y3=np.zeros(nodesperlayer*layers)
Z3=np.zeros(nodesperlayer*layers)
S3=np.zeros(nodesperlayer*layers)
index3=np.zeros((elementstot,4))
#Get new coordinates in 3d
for i in xrange(1,layers+1):
X3[i-1::layers]=X
Y3[i-1::layers]=Y
Z3[i-1::layers]=base+(i-1)*(surface-base)/(layers-1)
S3[i-1::layers]=S
if i<layers-1: #Build index3 with quads
ids=np.vstack((np.arange(i,nodestot-layers,layers),np.arange(i+1,nodestot-layers,layers),np.arange(i+layers+1,nodestot,layers),np.arange(i+layers,nodestot,layers))).T
index3[(i-1)*elementsperlayer:i*elementsperlayer,:]=ids
#Interpolation of data on specified points
data_interp=InterpFromMeshToMesh3d(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z,data,X3,Y3,Z3,np.nan)[0]
#build outputs
X=X3
Y=Y3
Z=Z3
S=S3
index=index3
return index,X,Y,Z,S,data_interp