<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 > 2)
+ continue;
+ end
+ planeSourcePoints(pointIndex,:) = [x,y];
+ offset = sum(planeSourcePoints(pointIndex,:).*planeNormal) - interceptDotNormal;
+ if(offset > 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) > 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 > 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 *, "ScalarRBFAndPlanar test"
call test_ibInterp_computeScalarRBFAndPlane
+print *, "VectorNoSlipRBFCoefficients test"
+ call test_ibInterp_computeVectorNoSlipRBFCoefficients
+
print *, "tests finished, stopping..."
stop
@@ -341,63 +344,132 @@
print *, "function normal deriv:", (functionAtDestination(1)-functionAtDestination(2))/(2.0*epsilon)
print *, "deriv should be:", 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, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients)
+
+ functionValue(1) = sum(coefficients(:,1)*testFunction)
+ functionValue(2) = sum(coefficients(:,1)*testFunction)
+ functionValue(3) = sum(coefficients(:,1)*testFunction)
+
+print *, "function value:", functionValue
+print *, "function value dot normal:", sum(functionValue*unitVectors(15,:))
+print *, "function should be:", testFunction(15)
+
+ end subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients
!
-! subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &
-! sourcePoints, unitVectors, destinationPoint, &
-! 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, &
-! sourcePoints, unitVectors, destinationPoint, &
-! alpha, matrix, rhs)
-!
-! do i = 1, pointCount
-! matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
-! matrix(pointCount+1:pointCount+3,i) &
-! = 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, &
! sourcePoints, unitVectors, destinationPoint, &
! alpha, planeBasisVectors, coefficients, dminfo)
</font>
</pre>