Index: /issm/trunk-jpl/src/m/extrusion/DepthAverage.m =================================================================== --- /issm/trunk-jpl/src/m/extrusion/DepthAverage.m (revision 17724) +++ /issm/trunk-jpl/src/m/extrusion/DepthAverage.m (revision 17725) @@ -1,4 +1,4 @@ function vector_average=DepthAverage(md,vector) -%DEPTHAVERAGE - computes depth average of 3d vector, and return value on 2d mesh. +%DEPTHAVERAGE - computes depth average of 3d vector using the trapezoidal rule, and returns the value on 2d mesh. % % Usage: Index: /issm/trunk-jpl/src/m/extrusion/DepthAverage.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/DepthAverage.py (revision 17725) +++ /issm/trunk-jpl/src/m/extrusion/DepthAverage.py (revision 17725) @@ -0,0 +1,50 @@ +import numpy as npy +from project2d import project2d + +def DepthAverage(md,vector): + ''' + computes depth average of 3d vector using the trapezoidal rule, and returns + the value on the 2d mesh. + + Usage: + vector_average=DepthAverage(md,vector) + + Example: + vel_bar=DepthAverage(md,md.initialization.vel) + ''' + + #check that the model given in input is 3d + if md.mesh.domaintype().lower() != '3d': + raise TypeError('DepthAverage error message: the model given in input must be 3d') + + # coerce to array in case float is passed + if type(vector)!=npy.ndarray: + print 'coercing array' + vector=npy.array(value) + + vec2d=False + if vector.ndim==2: + vec2d=True + vector=vector.reshape(-1,) + + #nods data + if vector.shape[0]==md.mesh.numberofvertices: + vector_average=npy.zeros(md.mesh.numberofvertices2d) + for i in xrange(1,md.mesh.numberoflayers): + vector_average=vector_average+(project2d(md,vector,i)+project2d(md,vector,i+1))/2.*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i)) + vector_average=vector_average/project2d(md,md.geometry.thickness,1) + + #element data + elif vector.shape[0]==md.mesh.numberofelements: + vector_average=npy.zeros(md.mesh.numberofelements2d) + for i in xrange(1,md.mesh.numberoflayers): + vector_average=vector_average+project2d(md,vector,i)*(project2d(md,md.mesh.z,i+1)-project2d(md,md.mesh.z,i)) + vector_average=vector_average/project2d(md,md.geometry.thickness,1) + + else: + raise ValueError('vector size not supported yet'); + + if vec2d: + vector_average=vector_average.reshape(-1,1) + + return vector_average Index: /issm/trunk-jpl/src/m/extrusion/project2d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project2d.py (revision 17725) +++ /issm/trunk-jpl/src/m/extrusion/project2d.py (revision 17725) @@ -0,0 +1,53 @@ +import numpy as npy + +def project2d(md3d,value,layer): + ''' + returns the value of a field for a given layer of the mesh + + project 'value' vector taken at layer 'layer' from extruded 2d-3d + mesh onto 2d mesh used to do the extrusion. This routine is used to + compare values between a 2d-3d mesh at a certain layer, and the + equivalent value (if it exists), on the original 2d mesh. + + Note: for consistency with usage in ISSM, layers are indexed starting + from one. This routine accounts for the fact that python indexes from + zero. In other words, when using this function consider the bottom + layer as layer=1. + + Usage: + projection_value=project2d(md3d,value,layer) + + Example: + vel2=project2d(md3d,md3d.vel,2) + ''' + + if md3d.mesh.domaintype().lower() != '3d': + raise StandardError("model passed to project2d function should be 3D") + + if layer<1 or layer>md3d.mesh.numberoflayers: + raise ValueError("layer must be between 0 and %i" % md3d.mesh.numberoflayers) + + # coerce to array in case float is passed + if type(value)!=npy.ndarray: + print 'coercing array' + value=npy.array(value) + + vec2d=False + if value.ndim==2: + value=value.reshape(-1,) + vec2d=True + + if value.size==1: + projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d] + elif value.shape[0]==md3d.mesh.numberofvertices: + #print 'indices: ', (layer-1)*md3d.mesh.numberofvertices2d, layer*md3d.mesh.numberofvertices2d + projection_value=value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d] + elif value.shape[0]==md3d.mesh.numberofvertices+1: + projection_value=[value[(layer-1)*md3d.mesh.numberofvertices2d:layer*md3d.mesh.numberofvertices2d], value[-1]] + else: + projection_value=value[(layer-1)*md3d.mesh.numberofelements2d:layer*md3d.mesh.numberofelements2d] + + if vec2d: + projection_value=projection_value.reshape(-1,1) + + return projection_value Index: /issm/trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 17724) +++ /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 17725) @@ -1,5 +1,4 @@ import numpy from pairoptions import pairoptions -import MatlabFuncs as m def project3d(md,*args): @@ -16,10 +15,10 @@ 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the value of the 2d vector. - 'padding': default to 0 (value adopted by other 3d layers not being projected0 + 'padding': default to 0 (value adopted by other 3d layers not being projected - Egs: - extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN); - extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0); - extruded_vector=project3d(md,'vector',vector2d,'type','node'); + Examples: + extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN) + extruded_vector=project3d(md,'vector',vector2d,'type','element','padding',0) + extruded_vector=project3d(md,'vector',vector2d,'type','node') """ @@ -27,5 +26,5 @@ if not md: raise TypeError("bad usage") - if not m.strcmp(md.mesh.elementtype(),'Penta'): + if md.mesh.domaintype().lower() != '3d': raise TypeError("input model is not 3d") @@ -33,5 +32,5 @@ options = pairoptions(*args) vector2d = options.getfieldvalue('vector') #mandatory - type = options.getfieldvalue('type') #mandatory + vectype = options.getfieldvalue('type') #mandatory layer = options.getfieldvalue('layer',0) #optional (do all layers otherwise) paddingvalue = options.getfieldvalue('padding',0) #0 by default @@ -40,15 +39,15 @@ if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1: vector1d=True - vector2d=vector2d.reshape(numpy.size(vector2d),1) + vector2d=vector2d.reshape(-1,1) - if isinstance(vector2d,(bool,int,long,float)) or numpy.size(vector2d)==1: + if isinstance(vector2d,(bool,int,long,float)) or numpy.size(vector2d)==1: projected_vector=vector2d - elif m.strcmpi(type,'node'): + elif vectype.lower()=='node': #Initialize 3d vector - if numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d: + if vector2d.shape[0]==md.mesh.numberofvertices2d: projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1: + elif vector2d.shape[0]==md.mesh.numberofvertices2d+1: projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) projected_vector[-1,:]=vector2d[-1,:] @@ -64,10 +63,10 @@ projected_vector[((layer-1)*md.mesh.numberofvertices2d):(layer*md.mesh.numberofvertices2d),:]=vector2d - elif m.strcmpi(type,'element'): + elif vectype.lower()=='element': #Initialize 3d vector - if numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d: + if vector2d.shape[0]==md.mesh.numberofelements2d: projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements, numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) - elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1: + elif vector2d.shape[0]==md.mesh.numberofelements2d+1: projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype) projected_vector[-1,:]=vector2d[-1,:] @@ -87,6 +86,5 @@ if vector1d: - projected_vector=projected_vector.reshape((numpy.size(projected_vector))) + projected_vector=projected_vector.reshape(-1,) return projected_vector -