<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)<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>