<p><b>ringler@lanl.gov</b> 2012-11-02 11:32:10 -0600 (Fri, 02 Nov 2012)</p><p>code to compute generalized barycentric coordinates<br>
</p><hr noshade><pre><font color="gray">Added: trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/areaBasedonVertexCoords.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/areaBasedonVertexCoords.m                                (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/areaBasedonVertexCoords.m        2012-11-02 17:32:10 UTC (rev 2295)
@@ -0,0 +1,17 @@
+function [ area ] = areaBasedonVertexCoords( tri )
+%areaBasedonVertexCoords : given vertex coord, compute area
+%coords need to be listed in CCW direction to get positive area
+
+x1=tri(1,1);
+x2=tri(2,1);
+x3=tri(3,1);
+
+y1=tri(1,2);
+y2=tri(2,2);
+y3=tri(3,2);
+
+area = 0.5*((x1-x2)*(y1-y3)-(y1-y2)*(x1-x3));
+
+
+end
+

Added: trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/testFunction.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/testFunction.m                                (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/testFunction.m        2012-11-02 17:32:10 UTC (rev 2295)
@@ -0,0 +1,11 @@
+function [ eval ] = testFunction( x,y )
+%testFunction: a linear test function
+
+c0 = 1.0;
+mx = 2.0;
+my = 5.0;
+
+eval = c0 + mx*x + my*y;
+
+end
+

Added: trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/wachspressCoordinates.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/wachspressCoordinates.m                                (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/Wachspress/wachspressCoordinates.m        2012-11-02 17:32:10 UTC (rev 2295)
@@ -0,0 +1,122 @@
+clear all
+
+%Wachspress coordinates
+
+nY=501;
+nX=501;
+
+gammaField=NaN(nX,nY);
+
+nTri=3 % number of vertices on a triangle
+nDim=2 % number of spatial dimensions
+
+%define a convex polygon by specifying the (x,y) of the vertices
+%nVertices=4
+%x=transpose(([-1 1 1 -1]))
+%y=transpose(([-1 -1 1 1]))
+
+nVertices=6
+b = sqrt(3)/2;
+x=transpose(([-1/2 1/2 1 1/2 -1/2 -1]))
+y=transpose(([-b -b 0 b b 0]))
+
+tri = zeros(nVertices,nTri,nDim)
+
+%find the B's of the Wachspress coordinate
+for i=1:nVertices
+    iP1=i+1
+    iM1=i-1
+    if i==nVertices
+        iP1=1
+    end
+    if i==1
+        iM1=nVertices
+    end
+    
+    tri(1,1,1)=x(i);
+    tri(1,1,2)=y(i);
+    tri(1,2,1)=x(iP1);
+    tri(1,2,2)=y(iP1);
+    tri(1,3,1)=x(iM1);
+    tri(1,3,2)=y(iM1);
+    
+    work(:,:) = tri(1,:,:);
+    [ B(i,1) ] = areaBasedonVertexCoords( work )
+    
+end
+
+%evaluate (linear) function at each vertex point
+for i=1:nVertices
+    [ q(i,1) ] = testFunction( x(i,1), y(i,1) );  
+end
+
+for jP=1:nY
+    jP
+    for iP=1:nX
+        
+        yP=-b + b*(jP-1)/((nY-1)/2);
+        xP=-1.0+(iP-1)/((nX-1)/2);
+        
+        rFrac = abs(yP)/b;
+        xEdge = 1.0 - rFrac/2.0;
+        
+        if abs(xP)&lt;xEdge
+        
+        %find triangle area, tri(nTri,nVert,x)
+        for i=1:nVertices
+          iP1=i+1;
+          if i==nVertices
+              iP1=1;
+          end
+          tri(i,1,1)=xP;
+          tri(i,1,2)=yP;
+          tri(i,2,1)=x(i,1);
+          tri(i,2,2)=y(i,1);
+          tri(i,3,1)=x(iP1,1);
+          tri(i,3,2)=y(iP1,1);
+    
+          % find area of triangles (the A's of Wachspress coordinates)
+          work(:,:) = tri(i,:,:);
+          [ A(i,1) ] = areaBasedonVertexCoords( work );
+    
+        end
+
+        %find the unnormalized wachspress weights
+        for i=1:nVertices
+    
+         iM1=i-1;
+         if i==1
+           iM1=nVertices;
+         end
+   
+         AMask = A;
+         AMask(i)=1.0;
+         AMask(iM1)=1.0;
+    
+         w(i,1)=1.0;
+         for j=1:nVertices
+             w(i,1) = w(i,1)*AMask(j,1);
+         end
+    
+         w(i,1) = B(i,1)*w(i,1);
+        end
+
+        total = sum(w(:,1));
+        gamma(:,1) = w(:,1) / total;
+
+        reconstruct = 0.0;
+        for i=1:nVertices
+            reconstruct = reconstruct + gamma(i,1)*q(i,1);
+        end
+        
+        reconstruct;
+        [ eval ] = testFunction( xP, yP );
+        
+        gammaField(iP,jP)=gamma(1,1);
+        
+        end
+        
+    end
+end
+
+contourf(gammaField)
\ No newline at end of file

</font>
</pre>