<p><b>xylar@lanl.gov</b> 2010-05-24 13:22:28 -0600 (Mon, 24 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Got the free-slip reconstruction to work<br>
<br>
Updated teh matlab test code, which now produces the same results as the rbfInterp_*_tanNeu_* routines<br>
<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m 2010-05-24 18:41:08 UTC (rev 303)
+++ branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m 2010-05-24 19:22:28 UTC (rev 304)
@@ -3,7 +3,7 @@
pointCount = length(functionSamples);
-coeffs = matrix\rhs;
+coeffs = matrix'\rhs;
functionValue = zeros(dim,1);
Modified: branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m 2010-05-24 18:41:08 UTC (rev 303)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m 2010-05-24 19:22:28 UTC (rev 304)
@@ -1,46 +1,106 @@
-clear all;
+clear variables;
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;
+[X,Y,Z] = ndgrid(-0.5:0.5,-0.5:0.5,0);
+
+index1 = 6;
+index2 = 8;
+
+sourcePoints = [X(:), Y(:), Z(:); X(:), Y(:), Z(:); 0,0,-0.5;0,0,0.5];
+
+normals = [ones(size(X(:))), zeros(size(Y(:))), zeros(size(Z(:))); zeros(size(X(:))), ones(size(Y(:))), zeros(size(Z(:))); 0,0,1;0,0,1];
+
normalIndices = zeros(size(sourcePoints,1),1);
isTangent = false(size(normalIndices));
-normalIndices(5) = 4;
-normalIndices(6) = 4;
-isTangent(5) = true;
-isTangent(6) = true;
+normalIndices(index1) = index1-length(X(:));
+normalIndices(index2) = index2-length(X(:));
+isTangent(index1) = true;
+isTangent(index2) = true;
-quiver3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), normals(:,1), normals(:,2), normals(:,3)); axis equal;
+figure(5);
+quiver3(sourcePoints(~isTangent,1), sourcePoints(~isTangent,2), sourcePoints(~isTangent,3),...
+ 0.8*normals(~isTangent,1), 0.8*normals(~isTangent,2), 0.8*normals(~isTangent,3), 0, 'b'); axis equal;
+hold on;
+quiver3(sourcePoints(isTangent,1), sourcePoints(isTangent,2), sourcePoints(isTangent,3),...
+ 0.8*normals(isTangent,1), 0.8*normals(isTangent,2), 0.8*normals(isTangent,3), 0, 'r'); axis equal;
+hold off;
-freeSlipFunction = zeros(size(sourcePoints));
-freeSlipFunction(:,2) = 1 + sourcePoints(:,1).^2;
-freeSlipFunction(5,2) = 2*sourcePoints(5,1);
-freeSlipFunction(6,2) = 0;
+freeSlipFunction = ones(size(sourcePoints));
+freeSlipFunction(:,2) = 1 + sourcePoints(:,1).*sourcePoints(:,2);
+freeSlipFunction(index1,2) = sourcePoints(index1,2);
+freeSlipFunction(index2,2) = sourcePoints(index2,2);
+
freeSlipSamples = sum(freeSlipFunction.*normals,2);
+matrix = constructRBFMatrix(sourcePoints, normals, normalIndices, isTangent, 3);
+cs = matrix\[freeSlipSamples; zeros(3,1)];
+
+% reconstructs the samples?
+
+% matrixTest = zeros(size(matrix));
+%
+% for(i=1:size(sourcePoints,1))
+% destPoint = sourcePoints(i,:);
+% if(isTangent(i))
+% normalVector = normals(normalIndices(i),:)
+% rhs = constructRBF_DerivRHS(destPoint, sourcePoints, normals, 2, normalVector);
+% else
+% rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 2);
+% end
+% for(j=1:size(matrix,2))
+% matrixTest(i,j) = sum(rhs(j,:).*normals(i,:));
+% end
+% end
+%
+% return;
+
[X,Y] = ndgrid(linspace(-0.5,0.5,21),linspace(-0.5,0.5,21));
functionValue = zeros(size(X,1),size(X,2),3);
+%nullSpace = null(matrix)
-matrix = constructRBFMatrix(sourcePoints, normals, normalIndices, isTangent, 3);
+%return;
-null(matrix)
+for(yIndex = 1:size(X,2))
+ for(xIndex = 1:size(X,1))
+ destPoint = [X(xIndex,yIndex),Y(xIndex,yIndex),0];
+ rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 3);
+ %functionValue(xIndex,yIndex,:) = rhs'*nullSpace;
+ functionValue(xIndex,yIndex,:) = rhs'*cs; %evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 2);
+ end
+end
-return
+figure(1);
+surface(X,Y,functionValue(:,:,2));
+surface(X,Y,1+X.*Y);
for(yIndex = 1:size(X,2))
for(xIndex = 1:size(X,1))
destPoint = [X(xIndex,yIndex),Y(xIndex,yIndex),0];
- rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 3);
- functionValue(xIndex,yIndex,:) = evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 3);
+ rhs = constructRBF_DerivRHS(destPoint, sourcePoints, normals, 3, [1,0,0]);
+ functionValue(xIndex,yIndex,:) = rhs'*cs; %evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 2);
end
end
+figure(2);
surface(X,Y,functionValue(:,:,2));
-surface(X,Y,1+X.^2);
+surface(X,Y,Y);
+
+epsilon = 1e-6;
+testIndex = 8;
+destPoint = sourcePoints(testIndex,:) - epsilon*normals(normalIndices(testIndex),:);
+rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 3);
+functionValue1 = evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 3);
+
+destPoint = sourcePoints(testIndex,:) + epsilon*normals(normalIndices(testIndex),:);
+rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 3);
+functionValue2 = evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 3);
+
+fnNormalDeriv = (functionValue2-functionValue1)'/(2*epsilon)
+
+tanFnNormalDeriv = sum(fnNormalDeriv.*normals(testIndex,:))
+
+freeSlipSamples(testIndex)
\ No newline at end of file
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-24 18:41:08 UTC (rev 303)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-24 19:22:28 UTC (rev 304)
@@ -1191,224 +1191,223 @@
end subroutine rbfInterp_func_3DPlane_vec_const_dir_compCoeffs
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! ! Purpose: Compute interpolation coefficients in 3D that can be used to
-! ! interpolate a number of vector functions at a given locations. This is useful
-! ! if the interpolation location does not change with time, or if several
-! ! fields are to be interpolated at a given time step. (If both the vector fields
-! ! and the interpolation locations vary with time, there is no clear advantage in
-! ! using either this method or the method for 2D interpoaltion above; for simplicity
-! ! and because we foresee more uses for the method of this subroutine, we have not
-! ! implemented a 3D version of the fixed field, variable interpolation location method
-! ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
-! ! Neumann tangential boundary conditions (such as free slip). The interpolation is
-! ! performed with basis functions that are RBFs plus a constant.
-! ! Input:
-! ! pointCount - the number of "source" points and functionValues supplied
-! ! sourcePoints - the location of the "source" points in the 3D space where the values of
-! ! the function are known
-! ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
-! ! tangent to the interface where the boundary condition will be applied. A Neumann
-! ! boundary condition will be applied at these points in these directions.
-! ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
-! ! gives the normal vector at the same sourcePoint. This information is needed to compute
-! ! the Neumann boundary condition at this point.
-! ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
-! ! is performed by supplying the value of the vector function dotted into each of these unit
-! ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
-! ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
-! ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
-! ! condition and the Neumann tangential boundary conditions at these points.
-! ! destinationPoint - the point where the interpolation will be performed
-! ! alpha - a constant that give the characteristic length scale of the RBFs,
-! ! should be on the order of the distance between points
-! ! Output:
-! ! coefficients - the coefficients used to interpolate a function with Dirichlet
-! ! boundary conditions to the specified destinationPoint
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
-! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
-! alpha, coefficients)
-!
-! integer, intent(in) :: pointCount
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
-! logical, dimension(pointCount), intent(in) :: isTangentToInterface
-! integer, dimension(pointCount), intent(in) :: normalVectorIndex
-! 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
-!
-! integer :: i, j
-! integer :: matrixSize
-!
-! real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
-! real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
-! integer, dimension(:), pointer :: pivotIndices
-!
-! matrixSize = pointCount+3 ! extra space for constant vector
-!
-! allocate(matrix(matrixSize,matrixSize))
-! allocate(matrixCopy(matrixSize,matrixSize))
-! allocate(rhs(matrixSize,3))
-! allocate(coeffs(matrixSize,3))
-! allocate(pivotIndices(matrixSize))
-!
-! matrix = 0.0
-! rhs = 0.0
-! coeffs = 0.0
-!
-! call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
-! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
-! alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-!
-! do i = 1, pointCount
-! if(.not. isTangentToInterface(i)) then
-! matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
-! end if
-! 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
-! matrixCopy = matrix
-! call LEGS(matrixCopy, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
-! end do
-! coefficients = coeffs(1:pointCount,:)
-!
-! deallocate(matrix)
-! deallocate(matrixCopy)
-! deallocate(rhs)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs
-!
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! ! Purpose: Compute interpolation coefficients in 3D that can be used to
-! ! interpolate a number of vector functions at a given locations. This is useful
-! ! if the interpolation location does not change with time, or if several
-! ! fields are to be interpolated at a given time step. (If both the vector fields
-! ! and the interpolation locations vary with time, there is no clear advantage in
-! ! using either this method or the method for 2D interpoaltion above; for simplicity
-! ! and because we foresee more uses for the method of this subroutine, we have not
-! ! implemented a 3D version of the fixed field, variable interpolation location method
-! ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
-! ! Neumann tangential boundary conditions (such as free slip). The interpolation is
-! ! performed with basis functions that are RBFs plus a constant.
-! ! Input:
-! ! pointCount - the number of "source" points and functionValues supplied
-! ! sourcePoints - the location of the "source" points in the 3D space where the values of
-! ! the function are known. The sourcePoints are projected into the plane given by
-! ! planeBasisVectors
-! ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
-! ! tangent to the interface where the boundary condition will be applied. A Neumann
-! ! boundary condition will be applied at these points in these directions.
-! ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
-! ! gives the normal vector at the same sourcePoint. This information is needed to compute
-! ! the Neumann boundary condition at this point.
-! ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
-! ! is performed by supplying the value of the vector function dotted into each of these unit
-! ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
-! ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
-! ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
-! ! condition and the Neumann tangential boundary conditions at these points. The unitVectors
-! ! are projected into the plane given by planeBasisVectors
-! ! destinationPoint - the point where the interpolation will be performed. The destinationPoint
-! ! is projected into the plane given by planeBasisVectors
-! ! alpha - a constant that give the characteristic length scale of the RBFs,
-! ! should be on the order of the distance between points
-! ! planeBasisVectors - the basis fectors for the plane where interpolation is performed.
-! ! All points are projected into this plane.
-! ! Output:
-! ! coefficients - the coefficients used to interpolate a function with Dirichlet
-! ! boundary conditions to the specified destinationPoint
-! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs(&
-! pointCount, sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, &
-! destinationPoint, alpha, planeBasisVectors, coefficients)
-!
-! integer, intent(in) :: pointCount
-! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
-! logical, dimension(pointCount), intent(in) :: isTangentToInterface
-! integer, dimension(pointCount), intent(in) :: normalVectorIndex
-! 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(2,3) :: planeBasisVectors
-! real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-!
-! integer :: i, j
-! integer :: matrixSize
-!
-! real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
-! real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
-! real(kind=RKIND), dimension(2) :: planarDestinationPoint
-!
-! real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
-! real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
-! integer, dimension(:), pointer :: pivotIndices
-!
-! matrixSize = pointCount+2 ! space for constant vector in plane
-!
-! allocate(matrix(matrixSize,matrixSize))
-! allocate(matrixCopy(matrixSize,matrixSize))
-! allocate(rhs(matrixSize,2))
-! allocate(coeffs(matrixSize,2))
-! allocate(pivotIndices(matrixSize))
-!
-! matrix = 0.0
-! rhs = 0.0
-! coeffs = 0.0
-!
-! do i = 1, pointCount
-! planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
-! planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
-! planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
-! planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
-! end do
-! planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
-! planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
-! call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
-! planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
-! planarDestinationPoint, alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-!
-! do i = 1, pointCount
-! if(.not. isTangentToInterface(i)) then
-! matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
-! matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
-! end if
-! matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
-! end do
-! do i = 1,2
-! rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
-! end do
-!
-! ! solve each linear system
-! matrixCopy = matrix
-! call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
-! call LEGS(matrixCopy, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
-!
-! coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
-! coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
-! coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
-! coefficients = coeffs(1:pointCount,:)
-!
-! deallocate(matrix)
-! deallocate(matrixCopy)
-! deallocate(rhs)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of vector functions at a given locations. This is useful
+ ! if the interpolation location does not change with time, or if several
+ ! fields are to be interpolated at a given time step. (If both the vector fields
+ ! and the interpolation locations vary with time, there is no clear advantage in
+ ! using either this method or the method for 2D interpoaltion above; for simplicity
+ ! and because we foresee more uses for the method of this subroutine, we have not
+ ! implemented a 3D version of the fixed field, variable interpolation location method
+ ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
+ ! Neumann tangential boundary conditions (such as free slip). The interpolation is
+ ! performed with basis functions that are RBFs plus a constant.
+ ! Input:
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! sourcePoints - the location of the "source" points in the 3D space where the values of
+ ! the function are known
+ ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
+ ! tangent to the interface where the boundary condition will be applied. A Neumann
+ ! boundary condition will be applied at these points in these directions.
+ ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
+ ! gives the normal vector at the same sourcePoint. This information is needed to compute
+ ! the Neumann boundary condition at this point.
+ ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
+ ! is performed by supplying the value of the vector function dotted into each of these unit
+ ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
+ ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
+ ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
+ ! condition and the Neumann tangential boundary conditions at these points.
+ ! destinationPoint - the point where the interpolation will be performed
+ ! alpha - a constant that give the characteristic length scale of the RBFs,
+ ! should be on the order of the distance between points
+ ! Output:
+ ! coefficients - the coefficients used to interpolate a function with Dirichlet
+ ! boundary conditions to the specified destinationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, coefficients)
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer, dimension(pointCount), intent(in) :: normalVectorIndex
+ 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
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
+ real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+3 ! extra space for constant vector
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(matrixCopy(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,3))
+ allocate(coeffs(matrixSize,3))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0
+ rhs = 0.0
+ coeffs = 0.0
+
+ call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+
+ do i = 1, pointCount
+ matrix(pointCount+1:pointCount+3,i) = unitVectors(i,:)
+ if(.not. isTangentToInterface(i)) then
+ matrix(i,pointCount+1:pointCount+3) = matrix(pointCount+1:pointCount+3,i)
+ end if
+ 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
+ matrixCopy = matrix
+ call LEGS(matrixCopy, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ end do
+ coefficients = coeffs(1:pointCount,:)
+
+ deallocate(matrix)
+ deallocate(matrixCopy)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of vector functions at a given locations. This is useful
+ ! if the interpolation location does not change with time, or if several
+ ! fields are to be interpolated at a given time step. (If both the vector fields
+ ! and the interpolation locations vary with time, there is no clear advantage in
+ ! using either this method or the method for 2D interpoaltion above; for simplicity
+ ! and because we foresee more uses for the method of this subroutine, we have not
+ ! implemented a 3D version of the fixed field, variable interpolation location method
+ ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
+ ! Neumann tangential boundary conditions (such as free slip). The interpolation is
+ ! performed with basis functions that are RBFs plus a constant.
+ ! Input:
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! sourcePoints - the location of the "source" points in the 3D space where the values of
+ ! the function are known. The sourcePoints are projected into the plane given by
+ ! planeBasisVectors
+ ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
+ ! tangent to the interface where the boundary condition will be applied. A Neumann
+ ! boundary condition will be applied at these points in these directions.
+ ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
+ ! gives the normal vector at the same sourcePoint. This information is needed to compute
+ ! the Neumann boundary condition at this point.
+ ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
+ ! is performed by supplying the value of the vector function dotted into each of these unit
+ ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
+ ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
+ ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
+ ! condition and the Neumann tangential boundary conditions at these points. The unitVectors
+ ! are projected into the plane given by planeBasisVectors
+ ! destinationPoint - the point where the interpolation will be performed. The destinationPoint
+ ! is projected into the plane given by planeBasisVectors
+ ! alpha - a constant that give the characteristic length scale of the RBFs,
+ ! should be on the order of the distance between points
+ ! planeBasisVectors - the basis fectors for the plane where interpolation is performed.
+ ! All points are projected into this plane.
+ ! Output:
+ ! coefficients - the coefficients used to interpolate a function with Dirichlet
+ ! boundary conditions to the specified destinationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs(&
+ pointCount, sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, &
+ destinationPoint, alpha, planeBasisVectors, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer, dimension(pointCount), intent(in) :: normalVectorIndex
+ 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(2,3), intent(in) :: planeBasisVectors
+ real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+ real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+ real(kind=RKIND), dimension(2) :: planarDestinationPoint
+
+ real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
+ real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+2 ! space for constant vector in plane
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(matrixCopy(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,2))
+ allocate(coeffs(matrixSize,2))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0
+ rhs = 0.0
+ coeffs = 0.0
+
+ do i = 1, pointCount
+ planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
+ planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
+ planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+ end do
+ planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
+ planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
+ call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
+ planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
+ planarDestinationPoint, alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+
+ do i = 1, pointCount
+ matrix(pointCount+1,i) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+ matrix(pointCount+2,i) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+ if(.not. isTangentToInterface(i)) then
+ matrix(i,pointCount+1:pointCount+2) = matrix(pointCount+1:pointCount+2,i)
+ end if
+ end do
+ do i = 1,2
+ rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+ end do
+
+ ! solve each linear system
+ matrixCopy = matrix
+ call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+ call LEGS(matrixCopy, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+
+ coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
+ coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
+ coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
+
+ deallocate(matrix)
+ deallocate(matrixCopy)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
+
+
!!!!!!!!!!!!!!!!!!!!!
! private subroutines
!!!!!!!!!!!!!!!!!!!!!
@@ -1584,6 +1583,7 @@
normalVector = unitVectors(normalVectorIndex(j),:)
normalDotX = sum(normalVector * (sourcePoints(j,:)-sourcePoints(i,:)))
call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
matrix(i,j) = rbfDerivOverR*normalDotX*unitVectorDotProduct
end do
</font>
</pre>