<p><b>ringler@lanl.gov</b> 2012-11-04 13:27:29 -0700 (Sun, 04 Nov 2012)</p><p><br>
renamed test directory Wachspress to WachspressCoordinates<br>
<br>
incorporated changes resulting from Rob's comments<br>
<br>
added text regarding forward trajectory calculation<br>
</p><hr noshade><pre><font color="gray">Added: trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/areaBasedonVertexCoords.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/areaBasedonVertexCoords.m (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/areaBasedonVertexCoords.m 2012-11-04 20:27:29 UTC (rev 2296)
@@ -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/WachspressCoordinates/basisFunction.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/basisFunction.m (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/basisFunction.m 2012-11-04 20:27:29 UTC (rev 2296)
@@ -0,0 +1,108 @@
+clear all
+
+%Wachspress coordinates
+
+%compute the basis function for one vertex
+
+%the basis function is plotted at the end
+
+nY=501;
+nX=501;
+
+gammaField=NaN(nX,nY);
+
+nTri=3 % number of vertices on a triangle
+nDim=2 % number of spatial dimensions
+
+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
+
+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;
+
+ gammaField(iP,jP)=gamma(1,1);
+
+ end
+
+ end
+end
+
+contourf(gammaField)
\ No newline at end of file
Added: trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/interpolateWachspress.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/interpolateWachspress.m (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/interpolateWachspress.m 2012-11-04 20:27:29 UTC (rev 2296)
@@ -0,0 +1,128 @@
+clear all
+
+%interpolate with Wachspress coordinates
+%given a hexagon, compare reconstruction of test function to exact soln
+
+% when test function is linear in x and y
+% reconstruction should be at the level of round off error
+
+% the routine plot the error at the of the computation
+
+nX=50;
+nY=50;
+reconstruct=NaN(nX,nY);
+eval=NaN(nX,nY);
+diff=NaN(nX,nY);
+
+nTri=3 % number of vertices on a triangle
+nDim=2 % number of spatial dimensions
+
+%build polygon (here it is a perfect hexagon)
+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]))
+
+%the polygon is broken into triangles
+%the number of triangles is nVertices
+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
+
+%scan across the polygon
+%if the coordinate is inside polygon, then reconstruct function
+%compare the reconstruction to the evaluation of the test function
+for jP=1:nY
+ 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(iP,jP)=0.0;
+ for i=1:nVertices
+ reconstruct(iP,jP) = reconstruct(iP,jP) + gamma(i,1)*q(i,1);
+ end
+
+ [ eval(iP,jP) ] = testFunction( xP, yP );
+
+ end
+
+ end
+end
+
+diff = reconstruct - eval;
+contourf(diff)
+colorbar
\ No newline at end of file
Added: trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/testFunction.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/testFunction.m (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/WachspressCoordinates/testFunction.m 2012-11-04 20:27:29 UTC (rev 2296)
@@ -0,0 +1,11 @@
+function [ eval ] = testFunction( x,y )
+%testFunction: a linear test function
+
+c0 = 2.0;
+mx = 2.0;
+my = 5.0;
+
+eval = c0 + mx*x + my*y;
+
+end
+
Modified: trunk/documents/ocean/current_design_doc/determining_pre_images/pre_images.pdf
===================================================================
(Binary files differ)
Modified: trunk/documents/ocean/current_design_doc/determining_pre_images/pre_images.tex
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/pre_images.tex 2012-11-02 17:32:10 UTC (rev 2295)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/pre_images.tex 2012-11-04 20:27:29 UTC (rev 2296)
@@ -185,7 +185,7 @@
\end{equation}
The $\lambda$ functions can be thought of as non-negative, basis functions for each vertex. An example of one of these basis function is shown in Figure \ref{fig:Basis}. These basis functions have some very desirable properties for our application. First, $\lambda_k({\bf x}_{v,j})=1$ for $k=j$ and zero otherwise. Along an edge connecting two vertices, say $k$ and $k+1$, the interpolation is a linear function of $\lambda_k$ and $\lambda_{k+1}$ only. This results in $\lambda_k=0$ along all edges not connected to ${\bf x}_{v,k}$. The interpolation is exact for linear functions (reference). And finally, if we impose no-slip boundary conditions as ${\bf R}_k=0$ at boundary vertices, this zero-velocity condition holds exactly along all boundary edges. (Note: A special fix will be needed for those parts of the ocean domain that are only one-edge).
-MATLAB code to test Wachspress coordinates is here:
+MATLAB code to test Wachspress coordinates is here: \\ https://www.dropbox.com/s/oggd4iyekeqg3ng/WachspressCoordinates.zip
\begin{figure}
\center{\includegraphics[width=14cm]{./vector_interpolation.pdf}}
@@ -234,6 +234,7 @@
\State ${\bf x}_{D}={\bf x}_A$-dtSegment*${\bf u}_{mid}$
\State find iCellBase; ${\bf x}_{D}$ resides in iCellBase
\State kTest = $\min$(maxLevelCell(iCellBase)-kLevel,kTest)
+ \State Linearly interpolate ${\bf u}^n_v$ and ${\bf u}^{n+1}_v$ to ${\bf u}^n_v(t_D)$
\State ${\bf u}_D$=vectorInterpolate(${\bf x}_{D},{\bf u}_v$(iCellBase,$t_D$))
\State ${\bf x}_{mid}=1/2({\bf x}_{D}+{\bf x}_{A})$
\State find iCellBase; ${\bf x}_{mid}$ resides in iCellBase
@@ -252,8 +253,8 @@
\section{Design Solution: Forward Trajectories}
Date last modified: 2011/01/05 \\
Contributors: Todd Ringler and Rob Lowrie \\
+The proposed algorithm is identical to the backward trajectory Algorithm \ref{alg:back_trajectories} with two exceptions. First, the negative sign used in the backward trajectory to indicate backward-in-time is replace with a positive sign. This occurs on line 22 and 32. Second, the starting location and velocity is the set to the quadrature point and velocity at quadrature point instead of the vertex location and vertex velocity. Given the similarity between the forward and backward trajectory calculation, it is anticipated that both will be done with exactly the same source code.
-
%-----------------------------------------------------------------------
\chapter{Design and Implementation}
Modified: trunk/documents/ocean/current_design_doc/determining_pre_images/vector_interpolation.key
===================================================================
(Binary files differ)
</font>
</pre>