Changeset 13482

Timestamp:

09/28/12 16:04:47 (12 years ago)

Author:

jschierm

Message:

NEW: New python functions in support of test514.py (which doesn't yet work).

Location:

issm/trunk-jpl/src/m

Files:

• mesh/ComputeHessian.m (modified) (2 diffs)
• mesh/ComputeHessian.py (added)
• mesh/ComputeMetric.m (modified) (4 diffs)
• mesh/ComputeMetric.py (added)
• mesh/GetNodalFunctionsCoeff.py (modified) (2 diffs)
• mesh/bamg.py (modified) (2 diffs)
• miscellaneous/MatlabFuncs.py (modified) (1 diff)

issm/trunk-jpl/src/m/mesh/ComputeHessian.m

@@ -20 +20 @@
         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
@@ -28 +28 @@
 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);
 

issm/trunk-jpl/src/m/mesh/ComputeMetric.m

@@ -9 +9 @@
 %      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);
@@ -18 +18 @@
 
 %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)
@@ -38 +38 @@
 
 %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)
@@ -54 +54 @@
                 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);

issm/trunk-jpl/src/m/mesh/GetNodalFunctionsCoeff.py

@@ -19 +19 @@
 
         #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
@@ -48 +48 @@
 
         #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

issm/trunk-jpl/src/m/mesh/bamg.py

@@ -262 +262 @@
 
         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...
@@ -310 +310 @@
 
         #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

issm/trunk-jpl/src/m/miscellaneous/MatlabFuncs.py

@@ -74 +74 @@
                 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
+