Source code for issm.averaging

import numpy as  np
from GetAreas import GetAreas
import MatlabFuncs as m
try:
	from scipy.sparse import csc_matrix
except ImportError:
	print "could not import scipy, no averaging capabilities enabled"

[docs]def averaging(md,data,iterations,layer=0): ''' AVERAGING - smooths the input over the mesh This routine takes a list over the elements or the nodes in input and return a list over the nodes. For each iterations it computes the average over each element (average of the vertices values) and then computes the average over each node by taking the average of the element around a node weighted by the elements volume For 3d mesh, a last argument can be added to specify the layer to be averaged on. Usage: smoothdata=averaging(md,data,iterations) smoothdata=averaging(md,data,iterations,layer) Examples: velsmoothed=averaging(md,md.initialization.vel,4) pressure=averaging(md,md.initialization.pressure,0) temperature=averaging(md,md.initialization.temperature,1,1) ''' if len(data)!=md.mesh.numberofelements and len(data)!=md.mesh.numberofvertices: raise StandardError('averaging error message: data not supported yet') if md.mesh.dimension()==3 and layer!=0: if layer<=0 or layer>md.mesh.numberoflayers: raise ValueError('layer should be between 1 and md.mesh.numberoflayers') else: layer=0 #initialization if layer==0: weights=np.zeros(md.mesh.numberofvertices,) data=data.flatten(1) else: weights=np.zeros(md.mesh.numberofvertices2d,) data=data[(layer-1)*md.mesh.numberofvertices2d+1:layer*md.mesh.numberofvertices2d,:] #load some variables (it is much faster if the variabes are loaded from md once for all) if layer==0: index=md.mesh.elements numberofnodes=md.mesh.numberofvertices numberofelements=md.mesh.numberofelements else: index=md.mesh.elements2d numberofnodes=md.mesh.numberofvertices2d numberofelements=md.mesh.numberofelements2d #build some variables if md.mesh.dimension()==3 and layer==0: rep=6 areas=GetAreas(index,md.mesh.x,md.mesh.y,md.mesh.z) elif md.mesh.dimension()==2: rep=3 areas=GetAreas(index,md.mesh.x,md.mesh.y) else: rep=3 areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d) index=index-1 # since python indexes starting from zero line=index.flatten(1) areas=np.vstack(areas).reshape(-1,) summation=1./rep*np.ones(rep,) linesize=rep*numberofelements #update weights that holds the volume of all the element holding the node i weights=csc_matrix( (np.tile(areas,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1)) #initialization if len(data)==numberofelements: average_node=csc_matrix( (np.tile(areas*data,(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1)) average_node=average_node/weights average_node = csc_matrix(average_node) else: average_node=csc_matrix(data.reshape(-1,1)) #loop over iteration for i in np.arange(1,iterations+1): average_el=np.asarray(np.dot(average_node.todense()[index].reshape(numberofelements,rep),np.vstack(summation))).reshape(-1,) average_node=csc_matrix( (np.tile(areas*average_el.reshape(-1),(rep,1)).reshape(-1,),(line,np.zeros(linesize,))), shape=(numberofnodes,1)) average_node=average_node/weights average_node=csc_matrix(average_node) #return output as a full matrix (C code do not like sparse matrices) average=np.asarray(average_node.todense()).reshape(-1,) return average