<p><b>xylar@lanl.gov</b> 2010-04-28 00:12:31 -0600 (Wed, 28 Apr 2010)</p><p>Branch Commit<br>
<br>
added matlab code to help with interpolation testing<br>
<br>
added a function for testing no-slip boundary condition in 3D<br>
  the no-slip code does not seem to be working properly :(<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/IBinterp/mpas/matlab/plane_test.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/plane_test.m                                (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/matlab/plane_test.m        2010-04-28 06:12:31 UTC (rev 218)
@@ -0,0 +1,49 @@
+sourcePoints = zeros(12,3);
+  
+planarBasisVectors = zeros(2,3);
+planarBasisVectors(1,:) = [0.816496580927726, -0.408248290463863, -0.408248290463863];
+planarBasisVectors(2,:) = [0, 0.707106781186548, -0.707106781186548];
+
+planeSourcePoints = zeros(12,2);
+planeNormal = [1,2];
+planeNormal = planeNormal/sqrt(sum(planeNormal.^2));
+
+planarIntercept = [1,1];
+interceptDotNormal = sum(planeNormal.*planarIntercept);
+
+
+testFunction = zeros(12,1);
+
+
+figure(1);
+pointIndex = 1;
+for(y = -1:2)
+    for(x = -1:2)
+        r = sqrt((x-0.5)^2 + (y-0.5)^2);
+        if(r &gt; 2)
+            continue;
+        end
+        planeSourcePoints(pointIndex,:) = [x,y];
+        offset = sum(planeSourcePoints(pointIndex,:).*planeNormal) - interceptDotNormal;
+        if(offset &gt; 0)
+            planeSourcePoints(pointIndex,:) = planeSourcePoints(pointIndex,:)...
+                - offset*planeNormal;
+            testFunction(pointIndex) = planeNormal(1)*2.0*planeSourcePoints(pointIndex,1) ...
+              *planeSourcePoints(pointIndex,2) + planeNormal(2)*planeSourcePoints(pointIndex,1)^2;
+plot(planeSourcePoints(pointIndex,1), planeSourcePoints(pointIndex,2), 'r.');
+hold on;
+        else
+            testFunction(pointIndex) = planeSourcePoints(pointIndex,1)^2*planeSourcePoints(pointIndex,2);
+plot(planeSourcePoints(pointIndex,1), planeSourcePoints(pointIndex,2), 'b.');            
+hold on;
+        end
+        sourcePoints(pointIndex,:) = planarBasisVectors(1,:)*planeSourcePoints(pointIndex,1)...
+            + planarBasisVectors(2,:)*planeSourcePoints(pointIndex,2);
+        pointIndex = pointIndex + 1;
+    end
+end
+hold off;
+axis equal;
+
+normal = planarBasisVectors(1,:)*planeNormal(1)...
+    + planarBasisVectors(2,:)*planeNormal(2);
\ No newline at end of file

Added: branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m                                (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m        2010-04-28 06:12:31 UTC (rev 218)
@@ -0,0 +1,85 @@
+planeNormal = [1,2,3];
+planeNormal = planeNormal./sqrt(sum(planeNormal.^2));
+
+planeTangents = zeros(2,3);
+planeTangents(1,:) = [1,0,0];
+planeTangents(1,:) = planeTangents(1,:) - sum(planeTangents(1,:).*planeNormal)*planeNormal;
+planeTangents(1,:) = planeTangents(1,:)/sqrt(sum(planeTangents(1,:).^2));
+planeTangents(2,:) = cross(planeNormal, planeTangents(1,:));
+
+planeIntercept = [0.5,0.5,0];
+
+normalDotIntercept = sum(planeNormal.*planeIntercept);
+
+edges = zeros(7^3*3,3);
+normals = zeros(size(edges));
+
+
+edgeIndex = 1;
+for(z = -3:3)
+    for(y = -3:3)
+        for(x = -3:3)
+            edges(edgeIndex,:) = [x+0.5,y,z];
+            normals(edgeIndex,:) = [1,0,0];
+            edges(edgeIndex+1,:) = [x,y+0.5,z];
+            normals(edgeIndex+1,:) = [0,1,0];
+            edges(edgeIndex+2,:) = [x,y,z+0.5];
+            normals(edgeIndex+2,:) = [0,0,1];
+            edgeIndex = edgeIndex + 3;
+        end
+    end
+end
+
+r = sqrt((edges(:,1)-0.0).^2 + (edges(:,2)-0.0).^2 + (edges(:,3)-0.0).^2);
+
+% edgeCumSum = cumsum(histc(r,0:0.1:2));
+% 
+% plot(0:0.1:2,edgeCumSum);
+
+indices = [];
+for(edgeIndex = 1:size(edges,1))
+    if(r(edgeIndex) &gt; 1.25)
+        continue;
+    end
+    indices(end+1) = edgeIndex;
+end
+
+sourcePoints = edges(indices,:);
+normals = normals(indices,:);
+
+offset = planeNormal(1)*sourcePoints(:,1)...
+    + planeNormal(2)*sourcePoints(:,2)...
+    + planeNormal(3)*sourcePoints(:,3) - normalDotIntercept;
+
+interfaceMask = offset &gt; 0;
+
+
+indices = find(interfaceMask);
+normalIndices = zeros(length(interfaceMask)+2*length(indices),1);
+isTangent = false(size(normalIndices));
+for(i = 1:length(indices))
+    index = indices(i);
+    sourcePoints(index,:) = sourcePoints(index,:)...
+        - offset(index)*planeNormal;
+    normals(index,:) = planeNormal;
+    sourcePoints(end+1,:) = sourcePoints(index,:);
+    normals(end+1,:) = planeTangents(1,:);
+    isTangent(size(normals,1)) = true;
+    normalIndices(size(normals,1)) = index;
+    sourcePoints(end+1,:) = sourcePoints(index,:);
+    normals(end+1,:) = planeTangents(2,:);
+    isTangent(size(normals,1)) = true;
+    normalIndices(size(normals,1)) = index;
+end
+
+plot3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), '.'); axis equal;
+
+noSlipFunction = [sourcePoints(:,1).^2.*sourcePoints(:,2),sourcePoints(:,3),sourcePoints(:,1)];
+freeSlipFunction = noSlipFunction;
+freeSlipFunction(interfaceMask,1) = planeNormal(1)*2*sourcePoints(interfaceMask,1).*sourcePoints(interfaceMask,2)...
+    + planeNormal(2)*sourcePoints(interfaceMask,1).^2;
+freeSlipFunction(interfaceMask,2) = planeNormal(3);
+freeSlipFunction(interfaceMask,3) = planeNormal(1);
+
+noSlipSamples = sum(noSlipFunction.*normals,2);
+freeSlipSamples = sum(freeSlipFunction.*normals,2);
\ No newline at end of file

Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F        2010-04-27 22:57:02 UTC (rev 217)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F        2010-04-28 06:12:31 UTC (rev 218)
@@ -32,6 +32,9 @@
 print *, &quot;ScalarRBFAndPlanar test&quot;
     call test_ibInterp_computeScalarRBFAndPlane
 
+print *, &quot;VectorNoSlipRBFCoefficients test&quot;
+    call test_ibInterp_computeVectorNoSlipRBFCoefficients
+
     print *, &quot;tests finished, stopping...&quot;
     stop
 
@@ -341,63 +344,132 @@
 print *, &quot;function normal deriv:&quot;, (functionAtDestination(1)-functionAtDestination(2))/(2.0*epsilon)
 print *, &quot;deriv should be:&quot;, testFunction(12)
   end subroutine test_ibInterp_computeScalarRBFAndPlane
+
+  subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients
+
+    integer, parameter :: pointCount = 46
+    real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+    real(kind=RKIND), dimension(pointCount,3) :: unitVectors
+    real(kind=RKIND), dimension(3) :: destinationPoint, functionValue
+    real(kind=RKIND) :: alpha
+    real(kind=RKIND), dimension(pointCount, 3) :: coefficients
+
+    real(kind=RKIND), dimension(pointCount) :: testFunction
+
+    alpha = 1.35
+
+    sourcePoints(1,:) = (/0.0, -0.500000000000000, -1.0/)
+    sourcePoints(2,:) = (/0.0, -1.0, -0.500000000000000/)
+    sourcePoints(3,:) = (/-0.500000000000000, 0.0, -1.0/)
+    sourcePoints(4,:) = (/-1.0, 0.0, -0.500000000000000/)
+    sourcePoints(5,:) = (/0.500000000000000, 0.0, -1.0/)
+    sourcePoints(6,:) = (/0.0, 0.500000000000000, -1.0/)
+    sourcePoints(7,:) = (/0.0, 0.0, -0.500000000000000/)
+    sourcePoints(8,:) = (/1.0, 0.0, -0.500000000000000/)
+    sourcePoints(9,:) = (/0.0, 1.0, -0.500000000000000/)
+    sourcePoints(10,:) = (/-0.500000000000000, -1.0, 0.0/)
+    sourcePoints(11,:) = (/-1.0, -0.500000000000000, 0.0/)
+    sourcePoints(12,:) = (/0.500000000000000, -1.0, 0.0/)
+    sourcePoints(13,:) = (/0.0, -0.500000000000000, 0.0/)
+    sourcePoints(14,:) = (/0.0, -1.0, 0.500000000000000/)
+    sourcePoints(15,:) = (/1.0, -0.500000000000000, 0.0/)
+    sourcePoints(16,:) = (/-0.500000000000000, 0.0, 0.0/)
+    sourcePoints(17,:) = (/-1.0, 0.500000000000000, 0.0/)
+    sourcePoints(18,:) = (/-1.0, 0.0, 0.500000000000000/)
+    sourcePoints(19,:) = (/0.500000000000000, 0.0, 0.0/)
+    sourcePoints(20,:) = (/0.0, 0.500000000000000, 0.0/)
+    sourcePoints(21,:) = (/0.0, 0.0, 0.500000000000000/)
+    sourcePoints(22,:) = (/0.964285714285714, 0.428571428571429, -0.107142857142857/)
+    sourcePoints(23,:) = (/0.928571428571429, -0.142857142857143, 0.285714285714286/)
+    sourcePoints(24,:) = (/-0.500000000000000, 1.0, 0.0/)
+    sourcePoints(25,:) = (/0.428571428571429, 0.857142857142857, -0.214285714285714/)
+    sourcePoints(26,:) = (/-0.142857142857143, 0.714285714285714, 0.0714285714285714/)
+    sourcePoints(27,:) = (/-0.0357142857142857, -0.571428571428571, 0.892857142857143/)
+    sourcePoints(28,:) = (/-0.571428571428571, -0.142857142857143, 0.785714285714286/)
+    sourcePoints(29,:) = (/0.357142857142857, -0.285714285714286, 0.571428571428571/)
+    sourcePoints(30,:) = (/-0.178571428571429, 0.142857142857143, 0.464285714285714/)
+    sourcePoints(31,:) = (/0.964285714285714, 0.428571428571429, -0.107142857142857/)
+    sourcePoints(32,:) = (/0.964285714285714, 0.428571428571429, -0.107142857142857/)
+    sourcePoints(33,:) = (/0.928571428571429, -0.142857142857143, 0.285714285714286/)
+    sourcePoints(34,:) = (/0.928571428571429, -0.142857142857143, 0.285714285714286/)
+    sourcePoints(35,:) = (/0.428571428571429, 0.857142857142857, -0.214285714285714/)
+    sourcePoints(36,:) = (/0.428571428571429, 0.857142857142857, -0.214285714285714/)
+    sourcePoints(37,:) = (/-0.142857142857143, 0.714285714285714, 0.0714285714285714/)
+    sourcePoints(38,:) = (/-0.142857142857143, 0.714285714285714, 0.0714285714285714/)
+    sourcePoints(39,:) = (/-0.0357142857142857, -0.571428571428571, 0.892857142857143/)
+    sourcePoints(40,:) = (/-0.0357142857142857, -0.571428571428571, 0.892857142857143/)
+    sourcePoints(41,:) = (/-0.571428571428571, -0.142857142857143, 0.785714285714286/)
+    sourcePoints(42,:) = (/-0.571428571428571, -0.142857142857143, 0.785714285714286/)
+    sourcePoints(43,:) = (/0.357142857142857, -0.285714285714286, 0.571428571428571/)
+    sourcePoints(44,:) = (/0.357142857142857, -0.285714285714286, 0.571428571428571/)
+    sourcePoints(45,:) = (/-0.178571428571429, 0.142857142857143, 0.464285714285714/)
+    sourcePoints(46,:) = (/-0.178571428571429, 0.142857142857143, 0.464285714285714/)
+
+    unitVectors(1,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(2,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(3,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(4,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(5,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(6,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(7,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(8,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(9,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(10,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(11,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(12,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(13,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(14,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(15,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(16,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(17,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(18,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(19,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(20,:) = (/0.0, 1.0, 0.0/)
+    unitVectors(21,:) = (/0.0, 0.0, 1.0/)
+    unitVectors(22,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(23,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(24,:) = (/1.0, 0.0, 0.0/)
+    unitVectors(25,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(26,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(27,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(28,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(29,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(30,:) = (/0.267261241912424, 0.534522483824849, 0.801783725737273/)
+    unitVectors(31,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(32,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(33,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(34,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(35,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(36,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(37,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(38,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(39,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(40,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(41,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(42,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(43,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(44,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+    unitVectors(45,:) = (/0.963624111659432, -0.148249863332220, -0.222374794998330/)
+    unitVectors(46,:) = (/1.38777878078145e-17, 0.832050294337844, -0.554700196225229/)
+
+    testFunction = (/-1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 1.0, 0.0, -0.250000000000000, 0.0, -0.250000000000000, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.822383489827826, 0.864313506097250, 0.250000000000000, 0.271157178383451, -0.0724644183610947, 0.448422287815161, -0.0506471741233457, 0.582052908771373, 0.106212968041363, 0.185460903190116, -0.624037720753383, -0.367545123997483, -0.277350098112615, 0.0881719303783468, -0.416025147168922, 0.0352255506168395, 0.138675049056307, -0.125126342375154, 0.762712769809690, -0.0343611199268557, 0.970725343394151, -0.199251274041264, 0.277350098112615, -0.0247308217712550, 0.485362671697076/)
+
+    destinationPoint = sourcePoints(15,:)
+
+    call ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &amp;
+      sourcePoints, unitVectors, destinationPoint, &amp;
+      alpha, coefficients)
+
+    functionValue(1) = sum(coefficients(:,1)*testFunction)
+    functionValue(2) = sum(coefficients(:,1)*testFunction)
+    functionValue(3) = sum(coefficients(:,1)*testFunction)
+
+print *, &quot;function value:&quot;, functionValue
+print *, &quot;function value dot normal:&quot;, sum(functionValue*unitVectors(15,:))
+print *, &quot;function should be:&quot;, testFunction(15)
+
+  end subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients
 !
-!  subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &amp;
-!    sourcePoints, unitVectors, destinationPoint, &amp;
-!    alpha, coefficients, dminfo)
-!
-!    integer, intent(in) :: pointCount
-!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
-!    real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
-!    real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
-!    real(kind=RKIND), intent(in) :: alpha
-!    real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-!    type (dm_info), intent(in) :: dminfo
-!
-!    integer :: i, j
-!    integer :: matrixSize
-!
-!    real(kind=RKIND), dimension(:,:), pointer :: matrix
-!    real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
-!    integer, dimension(:), pointer :: pivotIndices
-!
-!    matrixSize = pointCount+3 ! extra space for constant vector 
-!
-!    allocate(matrix(matrixSize,matrixSize))  
-!    allocate(rhs(matrixSize,3))  
-!    allocate(coeffs(matrixSize,3))  
-!    allocate(pivotIndices(matrixSize))  
-!
-!    matrix = 0.0
-!    rhs = 0.0
-!    coeffs = 0.0
-!
-!    call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 3, &amp;
-!      sourcePoints, unitVectors, destinationPoint, &amp;
-!      alpha, matrix, rhs)
-!
-!    do i = 1, pointCount
-!      matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
-!      matrix(pointCount+1:pointCount+3,i) &amp;
-!        = matrix(i,pointCount+1:pointCount+3)
-!    end do
-!    do i = 1, 3
-!      rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
-!    end do
-!
-!    ! solve each linear system
-!    do i = 1, 3
-!      call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
-!    end do
-!    coefficients = coeffs(1:pointCount,:)
-!
-!    deallocate(matrix)  
-!    deallocate(rhs)  
-!    deallocate(coeffs)  
-!    deallocate(pivotIndices) 
-!
-!  end subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients 
-!
 !  subroutine test_ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &amp;
 !    sourcePoints, unitVectors, destinationPoint, &amp;
 !    alpha, planeBasisVectors, coefficients, dminfo)

</font>
</pre>