import numpy as np
from issm.GetNodalFunctionsCoeff import GetNodalFunctionsCoeff
from issm.GetAreas import GetAreas
import issm.MatlabFuncs as m
[docs]def ComputeHessian(index,x,y,field,type):
"""
COMPUTEHESSIAN - compute hessian matrix from a field
Compute the hessian matrix of a given field
return the three components Hxx Hxy Hyy
for each element or each node
Usage:
hessian=ComputeHessian(index,x,y,field,type)
Example:
hessian=ComputeHessian(md.mesh.elements,md.mesh.x,md.mesh.y,md.inversion.vel_obs,'node')
"""
#some variables
numberofnodes=np.size(x)
numberofelements=np.size(index,axis=0)
#some checks
if np.size(field)!=numberofnodes and np.size(field)!=numberofelements:
raise TypeError("ComputeHessian error message: the given field size not supported yet")
if not m.strcmpi(type,'node') and not m.strcmpi(type,'element'):
raise TypeError("ComputeHessian error message: only 'node' or 'element' type supported yet")
#initialization
line=index.reshape(-1,order='F')
linesize=3*numberofelements
#get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma
[alpha,beta,dum]=GetNodalFunctionsCoeff(index,x,y)
areas=GetAreas(index,x,y)
#compute weights that hold the volume of all the element holding the node i
weights=m.sparse(line,np.ones((linesize,1)),np.tile(areas.reshape(-1,),(3,1)),numberofnodes,1)
#compute field on nodes if on elements
if np.size(field,axis=0)==numberofelements:
field=m.sparse(line,np.ones((linesize,1)),np.tile(areas*field,(3,1)),numberofnodes,1)/weights
#Compute gradient for each element
grad_elx=np.sum(field[index-1,0]*alpha,axis=1)
grad_ely=np.sum(field[index-1,0]*beta,axis=1)
#Compute gradient for each node (average of the elements around)
gradx=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_elx).reshape(-1,),(3,1)),numberofnodes,1)
grady=m.sparse(line,np.ones((linesize,1)),np.tile((areas*grad_ely).reshape(-1,),(3,1)),numberofnodes,1)
gradx=gradx/weights
grady=grady/weights
#Compute hessian for each element
hessian=np.vstack((np.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,),np.sum(grady[index-1,0]*beta,axis=1).reshape(-1,))).T
if m.strcmpi(type,'node'):
#Compute Hessian on the nodes (average of the elements around)
hessian=np.hstack((m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,0]).reshape(-1,),(3,1)),numberofnodes,1)/weights,
m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,1]).reshape(-1,),(3,1)),numberofnodes,1)/weights,
m.sparse(line,np.ones((linesize,1)),np.tile((areas*hessian[:,2]).reshape(-1,),(3,1)),numberofnodes,1)/weights ))
return hessian