<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} & n_{y,1} & n_{x,1} x_1 & n_{x,1} y_1 & n_{y,1} x_1 & n_{y,1} y_1 \\
+\begin{pmatrix} n_{x,1} & n_{x,1} x_1 & n_{x,1} y_1 & n_{y,1} & n_{y,1} x_1 & n_{y,1} y_1 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
- n_{x,k} & n_{y,k} & n_{x,k} x_k & n_{x,k} y_k & n_{y,k} x_k & n_{y,k} y_k
+ n_{x,k} & n_{x,k} x_k & n_{x,k} y_k & n_{y,k} & n_{y,k} x_k & 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>