<p><b>ringler@lanl.gov</b> 2012-11-04 15:16:13 -0700 (Sun, 04 Nov 2012)</p><p><br>
added matlab code for testing least squares reconstruction.<br>
<br>
adjusted least squares matrix entries in .tex source to match matlab code.<br>
</p><hr noshade><pre><font color="gray">Added: trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/leastSquares.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/leastSquares.m                                (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/leastSquares.m        2012-11-04 22:16:13 UTC (rev 2297)
@@ -0,0 +1,158 @@
+clear all
+
+%testing the least squares approach to vector reconstruction
+
+nVertices = 10;
+xVertex = zeros(nVertices,1);
+yVertex = zeros(nVertices,1);
+
+nEdges = 9;
+xEdge = zeros(nEdges,1);
+yEdge = zeros(nEdges,1);
+uVector = zeros(nEdges,2);
+normal = uVector;
+
+% everything between the ---- lines just builds the vertices and edges
+%{x,y}Vertex: position of vertices
+%{x,y}Edge: position of edges
+%normal: unit vector normal to edges
+%uVector: test vector evaluated at (xEdge,yEdge)
+%------------------------------------------------------------------------
+dx = 1.0;
+dtr = pi/180.0;
+delta = 120.0*dtr;
+
+xVertex(1)=0.0;
+yVertex(1)=0.0;
+
+angle = 0.0*dtr;
+iVertexBase=1;
+iCount = 1;
+iEdge = 0;
+for i=1:3
+    iCount=iCount+1;
+    xVertex(iCount) = xVertex(iVertexBase)+dx*cos(angle);
+    yVertex(iCount) = yVertex(iVertexBase)+dx*sin(angle);
+    iEdge=iEdge+1;
+    xEdge(iEdge) = 0.5*(xVertex(iVertexBase)+xVertex(iCount));
+    yEdge(iEdge) = 0.5*(yVertex(iVertexBase)+yVertex(iCount));
+    [uVector(iEdge,1), uVector(iEdge,2)] = ...
+                          testVector(xEdge(iEdge),yEdge(iEdge));
+    normal(iEdge,2) = +(xVertex(iCount) - xVertex(iVertexBase));
+    normal(iEdge,1) = -(yVertex(iCount) - yVertex(iVertexBase));
+    angle = angle + delta;
+end
+
+angle = 300.0*dtr;
+iVertexBase=2;
+for i=1:2
+    iCount=iCount+1;
+    xVertex(iCount) = xVertex(iVertexBase)+dx*cos(angle);
+    yVertex(iCount) = yVertex(iVertexBase)+dx*sin(angle);
+    iEdge=iEdge+1;
+    xEdge(iEdge) = 0.5*(xVertex(iVertexBase)+xVertex(iCount));
+    yEdge(iEdge) = 0.5*(yVertex(iVertexBase)+yVertex(iCount));
+    [uVector(iEdge,1), uVector(iEdge,2)] = ...
+                          testVector(xEdge(iEdge),yEdge(iEdge));
+    normal(iEdge,2) = +(xVertex(iCount) - xVertex(iVertexBase));
+    normal(iEdge,1) = -(yVertex(iCount) - yVertex(iVertexBase));
+    angle = angle + delta;
+end
+
+angle = 60.0*dtr;
+iVertexBase=3;
+for i=1:2
+    iCount=iCount+1;
+    xVertex(iCount) = xVertex(iVertexBase)+dx*cos(angle);
+    yVertex(iCount) = yVertex(iVertexBase)+dx*sin(angle);
+    iEdge=iEdge+1;
+    xEdge(iEdge) = 0.5*(xVertex(iVertexBase)+xVertex(iCount));
+    yEdge(iEdge) = 0.5*(yVertex(iVertexBase)+yVertex(iCount));
+    [uVector(iEdge,1), uVector(iEdge,2)] = ...
+                          testVector(xEdge(iEdge),yEdge(iEdge));
+    normal(iEdge,2) = +(xVertex(iCount) - xVertex(iVertexBase));
+    normal(iEdge,1) = -(yVertex(iCount) - yVertex(iVertexBase));
+    angle = angle + delta;
+end
+
+angle = 180.0*dtr;
+iVertexBase=4;
+for i=1:2
+    iCount=iCount+1;
+    xVertex(iCount) = xVertex(iVertexBase)+dx*cos(angle);
+    yVertex(iCount) = yVertex(iVertexBase)+dx*sin(angle);
+    iEdge=iEdge+1;
+    xEdge(iEdge) = 0.5*(xVertex(iVertexBase)+xVertex(iCount));
+    yEdge(iEdge) = 0.5*(yVertex(iVertexBase)+yVertex(iCount));
+    [uVector(iEdge,1), uVector(iEdge,2)] = ...
+                          testVector(xEdge(iEdge),yEdge(iEdge));
+    normal(iEdge,2) = +(xVertex(iCount) - xVertex(iVertexBase));
+    normal(iEdge,1) = -(yVertex(iCount) - yVertex(iVertexBase));
+    angle = angle + delta;
+end
+
+for iEdge=1:nEdges
+    mag = sqrt( normal(iEdge,1)^2 + normal(iEdge,2)^2);
+    normal(iEdge,:) = normal(iEdge,:) / mag;
+end
+%------------------------------------------------------------------------
+
+%set aside space for matrix and rhs
+M = zeros(nEdges,6);
+rhs = zeros(nEdges,1);
+
+%fill matrix and rhs
+%rhs is just the test vector dotted with the normal
+%Matrix(:,1) is n_x
+%Matrix(:,2) is n_x*x
+%Matrix(:,3) is n_x*y
+%Matrix(:,4) is n_x
+%Matrix(:,5) is n_y*x
+%Matrix(:,6) is n_y*y
+
+for iEdge=1:nEdges
+    rhs(iEdge) = normal(iEdge,1)*uVector(iEdge,1) + ...
+                 normal(iEdge,2)*uVector(iEdge,2);
+
+    M(iEdge,1) = normal(iEdge,1);
+    M(iEdge,2) = normal(iEdge,1)*xEdge(iEdge);
+    M(iEdge,3) = normal(iEdge,1)*yEdge(iEdge);
+    
+    M(iEdge,4) = normal(iEdge,2);
+    M(iEdge,5) = normal(iEdge,2)*xEdge(iEdge);
+    M(iEdge,6) = normal(iEdge,2)*yEdge(iEdge);
+    
+end
+
+%find psuedo inversion
+M = pinv(M);
+
+%solve!
+solution = M*rhs
+
+%put solution and uVector into single array so they are plotted to scale
+nTmp = nEdges+1;
+xTmp = zeros(nTmp,1);
+yTmp = zeros(nTmp,1);
+uTmp = zeros(nTmp,2);
+uTmp(1:nEdges,:) = uVector(:,:);
+xTmp(1:nEdges) = xEdge(:);
+yTmp(1:nEdges) = yEdge(:);
+uTmp(nTmp,1) = solution(1);
+uTmp(nTmp,2) = solution(4);
+
+scatter(xVertex,yVertex)
+axis([-3 3 -3 3])
+axis equal
+hold on
+scatter(xEdge,yEdge,'filled')
+quiver(xTmp,yTmp, uTmp(:,1),uTmp(:,2))
+%quiver(xEdge,yEdge,normal(:,1),normal(:,2))
+hold off
+
+    
+    
+
+    
+
+

Added: trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/testVector.m
===================================================================
--- trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/testVector.m                                (rev 0)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/leastSquaresReconstruction/testVector.m        2012-11-04 22:16:13 UTC (rev 2297)
@@ -0,0 +1,16 @@
+function [ u, v ] = testVector( x,y )
+%testVector: given x and y, produce a test vector
+
+u0 = 1.0;
+ux = 0.0;
+uy = -1.0;
+
+v0 = 1.0;
+vx = 1.0;
+vy = 0.0;
+
+u = u0 + ux*x + uy*y;
+v = v0 + vx*x + vy*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-04 20:27:29 UTC (rev 2296)
+++ trunk/documents/ocean/current_design_doc/determining_pre_images/pre_images.tex        2012-11-04 22:16:13 UTC (rev 2297)
@@ -145,11 +145,11 @@
 \end{figure} 
 
 \begin{equation}
-\begin{pmatrix} n_{x,1} &amp; n_{y,1} &amp; n_{x,1} x_1 &amp; n_{x,1} y_1 &amp; n_{y,1} x_1 &amp; n_{y,1} y_1  \\ 
+\begin{pmatrix} n_{x,1} &amp; n_{x,1} x_1 &amp; n_{x,1} y_1 &amp;  n_{y,1} &amp; n_{y,1} x_1 &amp; n_{y,1} y_1  \\ 
                           \vdots &amp; \vdots &amp; \vdots &amp; \vdots &amp; \vdots &amp; \vdots \\
-                          n_{x,k} &amp; n_{y,k} &amp; n_{x,k} x_k &amp; n_{x,k} y_k &amp; n_{y,k} x_k &amp; n_{y,k} y_k
+                          n_{x,k} &amp;  n_{x,k} x_k &amp; n_{x,k} y_k &amp; n_{y,k} &amp; n_{y,k} x_k &amp; n_{y,k} y_k
 \end{pmatrix}
-\begin{pmatrix} c_0 \\ d_0 \\ c_1 \\ c_2 \\ d_1 \\ d_2  \end{pmatrix} 
+\begin{pmatrix} c_0 \\ c_1 \\ c_2 \\  d_0 \\ d_1 \\ d_2  \end{pmatrix} 
 = 
 \begin{pmatrix} u_{n,1} \\ \vdots \\ u_{n,k} \end{pmatrix}
 \end{equation}
@@ -158,10 +158,9 @@
 \begin{equation}
 {\bf R}({\bf x}_v) = (c_0, d_0)
 \end{equation}
-MATLAB code that to test this least-squares reconstruction is here:
+MATLAB code to test this least-squares reconstruction is here: \\
+https://www.dropbox.com/s/69ql731ke1x8mwy/leastSquaresReconstruction.zip
 
-(Note: Having reconstructed full vectors at vertices, we can likely use this vector data to compute rate-of-strain with sufficient accuracy to be used in turbulence closure models.)
-
 \section{Design Solution: Vector Interpolation}
 Date last modified: 2011/01/05 \\
 Contributors: Todd Ringler and Pedro de Silva Peixoto \\

</font>
</pre>