<p><b>xylar@lanl.gov</b> 2010-05-13 06:33:31 -0600 (Thu, 13 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Work in progress on testing RBF interpolation for free-slip boundary condition<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m 2010-05-13 12:33:31 UTC (rev 268)
@@ -0,0 +1,71 @@
+function functionValue = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, functionSamples)
+
+% pointCount = size(sourcePoints,1);
+% matrixSize = pointCount+3;
+% matrix = zeros(matrixSize, matrixSize);
+% rhs = zeros(matrixSize,3);
+%
+% destPoint = sourcePoints(15,:);
+% alpha = 1.35;
+% for(j=1:pointCount)
+% for(i=1:pointCount)
+% rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
+% dotProduct = sum(normals(i,:).*normals(j,:));
+% matrix(i,j) = rbf*dotProduct;
+% end
+% rbf = 1/sqrt(1 + sum((destPoint-sourcePoints(j,:)).^2)/alpha^2);
+% rhs(j,:) = rbf*normals(j,:);
+% matrix(j,pointCount+(1:3)) = normals(j,:);
+% matrix(pointCount+(1:3),j) = normals(j,:);
+% end
+%
+% for(i=1:3)
+% rhs(pointCount+i,i) = 1;
+% end
+%
+% coeffs = matrix\rhs;
+%
+% functionValue(1) = sum(coeffs(1:pointCount,1).*noSlipSamples);
+% functionValue(2) = sum(coeffs(1:pointCount,2).*noSlipSamples);
+% functionValue(3) = sum(coeffs(1:pointCount,3).*noSlipSamples);
+
+pointCount = size(sourcePoints,1);
+matrixSize = pointCount+3;
+matrix = zeros(matrixSize, matrixSize);
+rhs = zeros(matrixSize,3);
+
+functionValue = zeros(1,3);
+
+alpha = 1.35;
+for(j=1:pointCount)
+ if(isTangent(j))
+ for(i=1:pointCount)
+ rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
+ normalVector = normals(normalIndices(j),:);
+ normalDotX = sum(normalVector.*(sourcePoints(j,:)-sourcePoints(i,:)));
+ rbfDerivOverR = -rbf^3;
+ dotProduct = sum(normals(i,:).*normals(j,:));
+ matrix(i,j) = rbfDerivOverR*dotProduct*normalDotX;
+ end
+ else
+ for(i=1:pointCount)
+ rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
+ dotProduct = sum(normals(i,:).*normals(j,:));
+ matrix(i,j) = rbf*dotProduct;
+ end
+ matrix(j,pointCount+(1:3)) = normals(j,:);
+ matrix(pointCount+(1:3),j) = normals(j,:);
+ end
+ rbf = 1/sqrt(1 + sum((destPoint-sourcePoints(j,:)).^2)/alpha^2);
+ rhs(j,:) = rbf*normals(j,:);
+end
+
+for(i=1:3)
+ rhs(pointCount+i,i) = 1;
+end
+
+coeffs = matrix\rhs;
+
+functionValue(1) = sum(coeffs(1:pointCount,1).*functionSamples);
+functionValue(2) = sum(coeffs(1:pointCount,2).*functionSamples);
+functionValue(3) = sum(coeffs(1:pointCount,3).*functionSamples);
\ No newline at end of file
Modified: branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m 2010-05-12 22:31:47 UTC (rev 267)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m 2010-05-13 12:33:31 UTC (rev 268)
@@ -72,43 +72,49 @@
normalIndices(size(normals,1)) = index;
end
-plot3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), '.'); axis equal;
+%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);
+freeSlipFunction(isTangent,1) = planeNormal(1)*2*sourcePoints(isTangent,1).*sourcePoints(isTangent,2)...
+ + planeNormal(2)*sourcePoints(isTangent,1).^2;
+freeSlipFunction(isTangent,2) = planeNormal(3);
+freeSlipFunction(isTangent,3) = planeNormal(1);
noSlipSamples = sum(noSlipFunction.*normals,2);
freeSlipSamples = sum(freeSlipFunction.*normals,2);
-pointCount = size(sourcePoints,1);
-matrixSize = pointCount+3;
-matrix = zeros(matrixSize, matrixSize);
-rhs = zeros(matrixSize,3);
-destPoint = sourcePoints(15,:);
-alpha = 1.35;
-for(j=1:pointCount)
- for(i=1:pointCount)
- rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
- dotProduct = sum(normals(i,:).*normals(j,:));
- matrix(i,j) = rbf*dotProduct;
+[X,Y] = ndgrid(linspace(0,1,21),linspace(0,1,21));
+
+functionValue = zeros(size(X,1),size(X,2),3);
+
+for(yIndex = 1:size(X,2))
+ for(xIndex = 1:size(X,1))
+ destPoint = [X(xIndex,yIndex),Y(xIndex,yIndex),0];
+ functionValue(xIndex,yIndex,:) = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, freeSlipSamples);
end
- rbf = 1/sqrt(1 + sum((destPoint-sourcePoints(j,:)).^2)/alpha^2);
- rhs(j,:) = rbf*normals(j,:);
- matrix(j,pointCount+(1:3)) = normals(j,:);
- matrix(pointCount+(1:3),j) = normals(j,:);
end
-for(i=1:3)
- rhs(pointCount+i,i) = 1;
-end
+close all;
+surface(X,Y,functionValue(:,:,1));
+surface(X,Y,X.^2.*Y);
-coeffs = matrix\rhs;
-
-functionValue(1) = sum(coeffs(1:pointCount,1).*noSlipSamples);
-functionValue(2) = sum(coeffs(1:pointCount,2).*noSlipSamples);
-functionValue(3) = sum(coeffs(1:pointCount,3).*noSlipSamples);
\ No newline at end of file
+% epsilon = 1e-6;
+%
+% functionValue = zeros(2,3);
+% testIndex = 31;
+%
+% destNormal = normals(normalIndices(testIndex),:);
+%
+% destPoint = sourcePoints(testIndex,:)+epsilon*destNormal;
+%
+% functionValue(1,:) = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, freeSlipSamples);
+%
+% destPoint = sourcePoints(testIndex,:)-epsilon*destNormal;
+%
+% functionValue(2,:) = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, freeSlipSamples);
+% functionNormalDeriv = (functionValue(1,:)-functionValue(2,:))/(2*epsilon);
+%
+% functionNormalDerivDotTan = sum(functionNormalDeriv.*normals(testIndex,:))
+% freeSlipSamples(testIndex)
Added: branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m 2010-05-13 12:33:31 UTC (rev 268)
@@ -0,0 +1,39 @@
+clear all;
+close all;
+
+sourcePoints = [-0.5, 0, 0; 0, -.5, 0; 0, .5, 0; .5, 0, 0; .5, 0, 0; .5, 0, 0; 0, 0, -.5; 0, 0, .5];
+%normals = [-1, 0, 0; 0, -1, 0; 0, 1, 0; 1, 0, 0; 0, 1, 0; 0, 0, 1; 0, 0, -1; 0, 0, 1];
+sourcePoints = sourcePoints([1:4,7:8],:);
+normals = zeros(size(sourcePoints));
+normals(:,2) = 1;
+
+normalIndices = zeros(size(sourcePoints,1),1);
+isTangent = false(size(normalIndices));
+normalIndices(5) = 4;
+normalIndices(6) = 4;
+isTangent(5) = true;
+isTangent(6) = true;
+
+%quiver3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), normals(:,1), normals(:,2), normals(:,3)); axis equal;
+
+freeSlipFunction = zeros(size(sourcePoints));
+freeSlipFunction(:,2) = 1 + sourcePoints(:,1).^2;
+freeSlipFunction(5,2) = 2*sourcePoints(5,1);
+freeSlipFunction(6,2) = 0;
+
+freeSlipSamples = sum(freeSlipFunction.*normals,2);
+
+
+[X,Y] = ndgrid(linspace(-0.5,0.5,21),linspace(-0.5,0.5,21));
+
+functionValue = zeros(size(X,1),size(X,2),3);
+
+for(yIndex = 1:size(X,2))
+ for(xIndex = 1:size(X,1))
+ destPoint = [X(xIndex,yIndex),Y(xIndex,yIndex),0];
+ functionValue(xIndex,yIndex,:) = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, freeSlipSamples);
+ end
+end
+
+surface(X,Y,functionValue(:,:,2));
+surface(X,Y,1+X.^2);
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-05-12 22:31:47 UTC (rev 267)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-05-13 12:33:31 UTC (rev 268)
@@ -18,6 +18,7 @@
use time_integration
use RBF_interpolation
use vector_reconstruction
+ use interp_tests
implicit none
@@ -33,6 +34,8 @@
call new_init_reconstruct(mesh)
+ call doTests(mesh)
+
call reconstruct(block % time_levs(1) % state, mesh)
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-12 22:31:47 UTC (rev 267)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-13 12:33:31 UTC (rev 268)
@@ -1579,7 +1579,7 @@
normalDotX = sum(normalVector * (sourcePoints(j,:)-sourcePoints(i,:)))
call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- matrix(i,j) = rbfDerivOverR*unitVectorDotProduct
+ matrix(i,j) = rbfDerivOverR*normalDotX*unitVectorDotProduct
end do
else
do i = 1, pointCount
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-05-12 22:31:47 UTC (rev 267)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-05-13 12:33:31 UTC (rev 268)
@@ -38,6 +38,9 @@
print *, "VectorDirichletPlanarRBFCoefficients test"
call test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
+print *, "VectorTangentNeumann test"
+ call test_VectorTangentNeumann
+
print *, "tests finished, stopping..."
stop
@@ -592,5 +595,168 @@
end subroutine test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
+ subroutine test_VectorTangentNeumann
+
+ integer, parameter :: pointCount = 46
+ real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,3) :: unitVectors
+ real(kind=RKIND), dimension(2,3) :: destinationPoint, functionValue
+ real(kind=RKIND) :: alpha, epsilon
+ real(kind=RKIND), dimension(pointCount, 3) :: coefficients
+
+ real(kind=RKIND), dimension(pointCount) :: testFunction
+ logical, dimension(pointCount) :: isTangentToInterface
+ integer, dimension(pointCount) :: normalVectorIndex
+
+ integer :: testIndex
+
+ alpha = 1.35
+ epsilon = 1e-6
+
+ 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.513512089453768, 0.518874521662771, 0.197500762285732, 0.518874521662771, 0.105522413421855, 0.518874521662771, -0.220343736838739, 0.518874521662771, -0.167127692138924, 0.518874521662771, 0.0319397343307530, 0.518874521662771, -0.165156727520413, 0.518874521662771, -0.175011550612971, 0.518874521662771/)
+
+ isTangentToInterface = .false.
+ isTangentToInterface(31:46) = .true.
+ normalVectorIndex = (/0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 22, 22, 23, 23, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30/)
+
+ testIndex = 22
+ destinationPoint(1,:) = sourcePoints(testIndex,:)
+
+ call rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint(1,:), &
+ alpha, coefficients)
+
+ functionValue(1,1) = sum(coefficients(:,1)*testFunction)
+ functionValue(1,2) = sum(coefficients(:,2)*testFunction)
+ functionValue(1,3) = sum(coefficients(:,3)*testFunction)
+
+print *, "function value:", functionValue(1,:)
+print *, "function should be about:", sourcePoints(testIndex,1)**2*sourcePoints(testIndex,2), sourcePoints(testIndex,3), sourcePoints(testIndex,1)
+print *, "function value dot normal:", sum(functionValue(1,:)*unitVectors(testIndex,:))
+print *, "function dot normal should be exactly:", testFunction(testIndex)
+
+ testIndex = 31
+ destinationPoint(1,:) = sourcePoints(testIndex,:)+epsilon*unitVectors(normalVectorIndex(testIndex),:)
+
+ call rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint(1,:), &
+ alpha, coefficients)
+
+ functionValue(1,1) = sum(coefficients(:,1)*testFunction)
+ functionValue(1,2) = sum(coefficients(:,2)*testFunction)
+ functionValue(1,3) = sum(coefficients(:,3)*testFunction)
+
+ destinationPoint(2,:) = sourcePoints(testIndex,:)-epsilon*unitVectors(normalVectorIndex(testIndex),:)
+
+ call rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint(2,:), &
+ alpha, coefficients)
+
+ functionValue(2,1) = sum(coefficients(:,1)*testFunction)
+ functionValue(2,2) = sum(coefficients(:,2)*testFunction)
+ functionValue(2,3) = sum(coefficients(:,3)*testFunction)
+
+print *, "function normal deriv:", (functionValue(2,:)-functionValue(1,:))/(2*epsilon)
+print *, "function normal deriv should be about:", unitVectors(normalVectorIndex(testIndex),1)*2*sourcePoints(testIndex,1)*sourcePoints(testIndex,2) + unitVectors(normalVectorIndex(testIndex),2)*sourcePoints(testIndex,1)**2, unitVectors(normalVectorIndex(testIndex),3), unitVectors(normalVectorIndex(testIndex),1)
+print *, "function normal deriv dot tangent:", sum(((functionValue(2,:)-functionValue(1,:))/(2*epsilon))*unitVectors(testIndex,:))
+print *, "function normal deriv dot tangent should be exactly:", testFunction(testIndex)
+
+ end subroutine test_VectorTangentNeumann
+
+
end module interp_tests
</font>
</pre>