Index: ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py =================================================================== --- ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 13481) +++ ../trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py (revision 13482) @@ -18,8 +18,8 @@ """ #make columns out of x and y - x=x.rshape(-1) - y=y.rshape(-1) + x=x.reshape(-1) + y=y.reshape(-1) #get nels and nods nels=numpy.size(index,axis=0) @@ -47,12 +47,12 @@ invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2)) #get alpha and beta - alpha=numpy.hstack((invdet*(y2-y3),invdet*(y3-y1),invdet*(y1-y2))) - beta =numpy.hstack((invdet*(x3-x2),invdet*(x1-x3),invdet*(x2-x1))) + alpha=numpy.hstack(((invdet*(y2-y3)).reshape(-1,1),(invdet*(y3-y1)).reshape(-1,1),(invdet*(y1-y2)).reshape(-1,1))) + beta =numpy.hstack(((invdet*(x3-x2)).reshape(-1,1),(invdet*(x1-x3)).reshape(-1,1),(invdet*(x2-x1)).reshape(-1,1))) #get gamma if requested gamma=numpy.zeros((nels,3)) - gamma=numpy.hstack((invdet*(x2*y3-x3*y2),invdet*(y1*x3-y3*x1),invdet*(x1*y2-x2*y1))) + gamma=numpy.hstack(((invdet*(x2*y3-x3*y2)).reshape(-1,1),(invdet*(y1*x3-y3*x1)).reshape(-1,1),(invdet*(x1*y2-x2*y1)).reshape(-1,1))) return alpha,beta,gamma Index: ../trunk-jpl/src/m/mesh/ComputeHessian.m =================================================================== --- ../trunk-jpl/src/m/mesh/ComputeHessian.m (revision 13481) +++ ../trunk-jpl/src/m/mesh/ComputeHessian.m (revision 13482) @@ -19,7 +19,7 @@ if length(field)~=numberofnodes & length(field)~=numberofelements, error('ComputeHessian error message: the given field size not supported yet'); end -if strcmpi(type,'node') & strcmpi(type,'element'), +if ~strcmpi(type,'node') & ~strcmpi(type,'element'), error('ComputeHessian error message: only ''node'' or ''element'' type supported yet'); end @@ -27,11 +27,11 @@ line=index(:); linesize=3*numberofelements; -%get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma +%get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma [alpha beta]=GetNodalFunctionsCoeff(index,x,y); areas=GetAreas(index,x,y); -%comput weights that holds the volume of all the element holding the node i +%compute weights that hold the volume of all the element holding the node i weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofnodes,1); %compute field on nodes if on elements Index: ../trunk-jpl/src/m/mesh/bamg.py =================================================================== --- ../trunk-jpl/src/m/mesh/bamg.py (revision 13481) +++ ../trunk-jpl/src/m/mesh/bamg.py (revision 13482) @@ -261,7 +261,7 @@ #bamg_geometry=processgeometry(bamg_geometry,options.getfieldvalue('tol',float(nan)),domain[0]) elif isinstance(md.private.bamg,dict) and 'geometry' in md.private.bamg: - bamg_geometry=bamggeom(md.private.bamg['geometry']) + bamg_geometry=bamggeom(md.private.bamg['geometry'].__dict__) else: #do nothing... pass @@ -309,7 +309,7 @@ #}}} #call Bamg - bamgmesh_out,bamggeom_out=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options) + [bamgmesh_out,bamggeom_out]=BamgMesher(bamg_mesh.__dict__,bamg_geometry.__dict__,bamg_options) # plug results onto model md.private.bamg=OrderedDict() Index: ../trunk-jpl/src/m/mesh/ComputeMetric.py =================================================================== --- ../trunk-jpl/src/m/mesh/ComputeMetric.py (revision 0) +++ ../trunk-jpl/src/m/mesh/ComputeMetric.py (revision 13482) @@ -0,0 +1,74 @@ +import numpy + +def ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos): + """ + COMPUTEMETRIC - compute metric from an Hessian + + Usage: + metric=ComputeMetric(hessian,scale,epsilon,hmin,hmax,pos) + pos is contains the positions where the metric is wished to be maximized (water?) + + Example: + metric=ComputeMetric(hessian,2/9,10^-1,100,10^5,find(md.nodeonwater) + """ + + #first, find the eigen values of each line of H=[hessian(i,1) hessian(i,2); hessian(i,2) hessian(i,3)] + a=hessian[:,0] + b=hessian[:,1] + d=hessian[:,2] + lambda1=0.5*((a+d)+numpy.sqrt(4.*b**2+(a-d)**2)) + lambda2=0.5*((a+d)-numpy.sqrt(4.*b**2+(a-d)**2)) + pos1=numpy.nonzero(lambda1==0)[0] + pos2=numpy.nonzero(lambda2==0)[0] + pos3=numpy.nonzero(numpy.logical_and(b==0,lambda1==lambda2))[0] + + #Modify the eigen values to control the shape of the elements + lambda1=numpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) + lambda2=numpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) + + #compute eigen vectors + norm1=numpy.sqrt(8.*b**2+2.*(d-a)**2+2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2)) + v1x=2.*b/norm1 + v1y=((d-a)+numpy.sqrt((a-d)**2+4.*b**2))/norm1 + norm2=numpy.sqrt(8.*b**2+2.*(d-a)**2-2.*(d-a)*numpy.sqrt((a-d)**2+4.*b**2)) + v2x=2.*b/norm2 + v2y=((d-a)-numpy.sqrt((a-d)**2+4.*b**2))/norm2 + + v1x[pos3]=1. + v1y[pos3]=0. + v2x[pos3]=0. + v2y[pos3]=1. + + #Compute new metric (for each node M=V*Lambda*V^-1) + metric=numpy.hstack((((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v2y*v1x-lambda2*v1y*v2x)).reshape(-1,1), \ + ((v1x*v2y-v1y*v2x)**(-1)*( lambda1*v1y*v2y-lambda2*v1y*v2y)).reshape(-1,1), \ + ((v1x*v2y-v1y*v2x)**(-1)*(-lambda1*v2x*v1y+lambda2*v1x*v2y)).reshape(-1,1))) + + #some corrections for 0 eigen values + metric[pos1,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos1),1)) + metric[pos2,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos2),1)) + + #take care of water elements + metric[pos ,:]=numpy.tile(numpy.array([[1./hmax**2,0.,1./hmax**2]]),(numpy.size(pos ),1)) + + #take care of NaNs if any (use Numpy eig in a loop) + pos=numpy.nonzero(numpy.isnan(metric))[0] + if pos: + print(" %i NaN found in the metric. Use Numpy routine..." % numpy.size(pos)) + for posi in pos: + H=numpy.array([[hessian[posi,0],hessian[posi,1]],[hessian[posi,1],hessian[posi,2]]]) + [v,u]=numpy.linalg.eig(H) + v=numpy.diag(v) + lambda1=v[0,0] + lambda2=v[1,1] + v[0,0]=numpy.minimum(numpy.maximum(numpy.abs(lambda1)*scale/epsilon,1./hmax**2),1./hmin**2) + v[1,1]=numpy.minimum(numpy.maximum(numpy.abs(lambda2)*scale/epsilon,1./hmax**2),1./hmin**2) + + metricTria=numpy.dot(numpy.dot(u,v),numpy.linalg.inv(u)) + metric[posi,:]=numpy.array([metricTria[0,0],metricTria[0,1],metricTria[1,1]]) + + if numpy.any(numpy.isnan(metric)): + raise RunTimeError("ComputeMetric error message: NaN in the metric despite our efforts...") + + return metric + Index: ../trunk-jpl/src/m/mesh/ComputeMetric.m =================================================================== --- ../trunk-jpl/src/m/mesh/ComputeMetric.m (revision 13481) +++ ../trunk-jpl/src/m/mesh/ComputeMetric.m (revision 13482) @@ -8,28 +8,28 @@ % Example: % metric=ComputeMetric(hessian,2/9,10^-1,100,10^5,find(md.nodeonwater) -%first, find the eigen values of eah line of H=[hessian(i,1) hessian(i,2); hessian(i,2) hessian(i,3)] +%first, find the eigen values of each line of H=[hessian(i,1) hessian(i,2); hessian(i,2) hessian(i,3)] a=hessian(:,1); b=hessian(:,2); d=hessian(:,3); -lambda1=0.5*((a+d)+sqrt(4*b.^2+(a-d).^2)); -lambda2=0.5*((a+d)-sqrt(4*b.^2+(a-d).^2)); +lambda1=0.5*((a+d)+sqrt(4.*b.^2+(a-d).^2)); +lambda2=0.5*((a+d)-sqrt(4.*b.^2+(a-d).^2)); pos1=find(lambda1==0); pos2=find(lambda2==0); pos3=find(b==0 & lambda1==lambda2); %Modify the eigen values to control the shape of the elements -lambda1=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2); -lambda2=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2); +lambda1=min(max(abs(lambda1)*scale/epsilon,1./hmax^2),1./hmin^2); +lambda2=min(max(abs(lambda2)*scale/epsilon,1./hmax^2),1./hmin^2); %compute eigen vectors -norm1=sqrt(8*b.^2+2*(d-a).^2+2*(d-a).*sqrt((a-d).^2+4*b.^2)); -v1x=2*b./norm1; -v1y=((d-a)+sqrt((a-d).^2+4*b.^2))./norm1; -norm2=sqrt(8*b.^2+2*(d-a).^2-2*(d-a).*sqrt((a-d).^2+4*b.^2)); -v2x=2*b./norm2; -v2y=((d-a)-sqrt((a-d).^2+4*b.^2))./norm2; +norm1=sqrt(8.*b.^2+2.*(d-a).^2+2.*(d-a).*sqrt((a-d).^2+4.*b.^2)); +v1x=2.*b./norm1; +v1y=((d-a)+sqrt((a-d).^2+4.*b.^2))./norm1; +norm2=sqrt(8.*b.^2+2.*(d-a).^2-2.*(d-a).*sqrt((a-d).^2+4.*b.^2)); +v2x=2.*b./norm2; +v2y=((d-a)-sqrt((a-d).^2+4.*b.^2))./norm2; -v1x(pos3)=1; v1y(pos3)=0; -v2x(pos3)=0; v2y(pos3)=1; +v1x(pos3)=1.; v1y(pos3)=0.; +v2x(pos3)=0.; v2y(pos3)=1.; %Compute new metric (for each node M=V*Lambda*V^-1) metric=full([(v1x.*v2y-v1y.*v2x).^(-1).*(lambda1.*v2y.*v1x-lambda2.*v1y.*v2x) ... @@ -37,11 +37,11 @@ (v1x.*v2y-v1y.*v2x).^(-1).*(-lambda1.*v2x.*v1y+lambda2.*v1x.*v2y)]); %some corrections for 0 eigen values -metric(pos1,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos1),1); -metric(pos2,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos2),1); +metric(pos1,:)=repmat([1./hmax^2 0. 1./hmax^2],length(pos1),1); +metric(pos2,:)=repmat([1./hmax^2 0. 1./hmax^2],length(pos2),1); %take care of water elements -metric(pos,:)=repmat([1/hmax^2 0 1/hmax^2],length(pos),1); +metric(pos,:)=repmat([1./hmax^2 0. 1./hmax^2],length(pos),1); %take care of NaNs if any (use Matlab eig in a loop) [pos posj]=find(isnan(metric)); clear posj; @@ -53,8 +53,8 @@ [u,v]=eig(full(H)); lambda1=v(1,1); lambda2=v(2,2); - v(1,1)=min(max(abs(lambda1)*scale/epsilon,1/hmax^2),1/hmin^2); - v(2,2)=min(max(abs(lambda2)*scale/epsilon,1/hmax^2),1/hmin^2); + v(1,1)=min(max(abs(lambda1)*scale/epsilon,1./hmax^2),1./hmin^2); + v(2,2)=min(max(abs(lambda2)*scale/epsilon,1./hmax^2),1./hmin^2); metricTria=u*v*u^(-1); metric(pos(i),:)=[metricTria(1,1) metricTria(1,2) metricTria(2,2)]; Index: ../trunk-jpl/src/m/mesh/ComputeHessian.py =================================================================== --- ../trunk-jpl/src/m/mesh/ComputeHessian.py (revision 0) +++ ../trunk-jpl/src/m/mesh/ComputeHessian.py (revision 13482) @@ -0,0 +1,66 @@ +import numpy +from GetNodalFunctionsCoeff import * +from GetAreas import * +from MatlabFuncs import * + +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=numpy.size(x) + numberofelements=numpy.size(index,axis=0) + + #some checks + if numpy.size(field)!=numberofnodes and numpy.size(field)!=numberofelements: + raise TypeError("ComputeHessian error message: the given field size not supported yet") + if not strcmpi(type,'node') and not 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=sparse(line,numpy.ones((linesize,1)),numpy.tile(areas.reshape(-1,1),(3,1)),numberofnodes,1) + + #compute field on nodes if on elements + if numpy.size(field,axis=0)==numberofelements: + field=sparse(line,numpy.ones((linesize,1)),numpy.tile(areas*field,(3,1)),numberofnodes,1)/weights + + #Compute gradient for each element + grad_elx=numpy.sum(field[index.astype(int)-1,0]*alpha,axis=1) + grad_ely=numpy.sum(field[index.astype(int)-1,0]*beta,axis=1) + + #Compute gradient for each node (average of the elements around) + gradx=sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*grad_elx).reshape(-1,1),(3,1)),numberofnodes,1) + grady=sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*grad_ely).reshape(-1,1),(3,1)),numberofnodes,1) + gradx=gradx/weights + grady=grady/weights + + #Compute hessian for each element + hessian=numpy.hstack((numpy.sum(gradx[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*beta,axis=1).reshape(-1,1))) + + if strcmpi(type,'node'): + #Compute Hessian on the nodes (average of the elements around) + hessian=numpy.hstack((sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,0]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \ + sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,1]).reshape(-1,1),(3,1)),numberofnodes,1)/weights, \ + sparse(line,numpy.ones((linesize,1)),numpy.tile((areas*hessian[:,2]).reshape(-1,1),(3,1)),numberofnodes,1)/weights )) + + return hessian + Index: ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py =================================================================== --- ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 13481) +++ ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 13482) @@ -73,3 +73,18 @@ else: raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape)) +def sparse(ivec,jvec,svec,m=0,n=0,nzmax=0): + import numpy + + if not m: + m=numpy.max(ivec) + if not n: + n=numpy.max(jvec) + + a=numpy.zeros((m,n)) + + for i,j,s in zip(ivec.reshape(-1,order='F'),jvec.reshape(-1,order='F'),svec.reshape(-1,order='F')): + a[i-1,j-1]+=s + + return a +