<p><b>xylar@lanl.gov</b> 2010-05-05 03:41:20 -0600 (Wed, 05 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Fully commented RBF interpolation functions, renamed to be as descriptive as possible<br>
<br>
Vector reconstruction has been replaced with the new code and seems to work nearly exactly as before.<br>
<br>
Free slip vector reconstruction still requires testing before the modules are ready to be integrated into the trunk.<br>
</p><hr noshade><pre><font color="gray">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-05 03:51:01 UTC (rev 248)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-05-05 09:41:20 UTC (rev 249)
@@ -16,7 +16,7 @@
use grid_types
use time_integration
- use IB_interpolation
+ use RBF_interpolation
use vector_reconstruction
implicit none
@@ -27,7 +27,7 @@
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
- call ibInterp_initialize(mesh)
+ call rbfInterp_initialize(mesh)
!call init_reconstruct(mesh)
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-05 03:51:01 UTC (rev 248)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-05 09:41:20 UTC (rev 249)
@@ -20,10 +20,10 @@
! for finding the location on the the RBF reconstruction of a function
! (e.g., a height field) that minimizes the distance to a point in 3D space
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- public :: rbfInterp_location_2D_scalar_constant_computeCoefficients, &
- rbfInterp_location_2D_scalar_linear_computeCoefficients, &
- rbfInterp_location_2D_scalar_constant_evaluateWithDerivs, &
- rbfInterp_location_2D_scalar_linear_evaluateWithDerivs
+ public :: rbfInterp_loc_2D_sca_const_compCoeffs, &
+ rbfInterp_loc_2D_sca_lin_compCoeffs, &
+ rbfInterp_loc_2D_sca_const_evalWithDerivs, &
+ rbfInterp_loc_2D_sca_lin_evalWithDerivs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Routines for computing scalar interpolaiton coefficients in 3D (coplanar points
@@ -32,17 +32,20 @@
! above, these coefficients are valid for computing the value of the
! interpolant at a fixe point in space using function values that may
! vary (e.g., in time) at each of the interpolation "source" points.
- ! If boundary points are provided, coefficients for imposing both Neumann
- ! and Dirichlet boundary conditions are provided
+ ! The last 3 routines can be used to compute coefficients for imposing both Neumann
+ ! and Dirichlet boundary conditions.
! Pseudocode for function reconstruction at the destinationPoint is as follows
! Dirichlet: functionAtDestination = sum(functionAtSources*dirichletCoefficients)
- ! Neumann: functionAtDestination &
- ! = sum(functionAtSources(.not.isInterface)*neumannCoefficients(.not.isInterface)) &
- ! + sum(functionNormalDerivAtSources(isInterface)*neumannCoefficients(isInterface))
+ ! Neumann: functionAtDestination = sum(functionOrDerivAtSources*neumannCoefficients)
+ ! where functionOrDerivAtSource = functionAtSources where .not.(isInterface)
+ ! = functionNormalDerivAtSources where isInterface
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- public :: rbfInterp_function_3D_scalar_constant_computeCoefficients, &
- rbfInterp_function_3DPlane_scalar_linear_computeCoefficients, &
- rbfInterp_function_3D_scalar_linear_computeCoefficients
+ public :: rbfInterp_func_3D_sca_const_dir_compCoeffs, &
+ rbfInterp_func_3DPlane_sca_lin_dir_compCoeffs, &
+ rbfInterp_func_3D_sca_lin_dir_compCoeffs, &
+ rbfInterp_func_3D_sca_const_dirNeu_compCoeffs, &
+ rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs, &
+ rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Routines for computing vector interpolaiton coefficients in 3D (coplanar points
@@ -83,10 +86,10 @@
! + sum_(j where isTangent) ((dFunction/dn dot UnitVector)_j*coefficients_j,i)
! for i = x,y,z
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- public :: rbfInterp_function_3D_vector_constant_dirichlet_computeCoefficients, &
- rbfInterp_function_3DPlane_vector_constant_dirichlet_computeCoefficients, &
- rbfInterp_function_3D_vector_constant_tangentNeumann_computeCoefficients, &
- rbfInterp_function_3DPlane_vector_constant_tangentNeumann_computeCoefficients
+ public :: rbfInterp_func_3D_vec_const_dir_compCoeffs, &
+ rbfInterp_func_3DPlane_vec_const_dir_compCoeffs, &
+ rbfInterp_func_3D_vec_const_tanNeu_compCoeffs, &
+ rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
contains
@@ -143,7 +146,7 @@
radialUnitVectors(1,iCell) = xCell(iCell)
radialUnitVectors(2,iCell) = yCell(iCell)
radialUnitVectors(3,iCell) = zCell(iCell)
- call unit_vector_in_R3(radialUnitVectors(:,iCell))
+ call unit_vec_in_R3(radialUnitVectors(:,iCell))
end do
do iEdge = 1,nEdges
@@ -153,7 +156,7 @@
edgeNormalVectors(1,iEdge) = xCell(cell2) - xCell(cell1)
edgeNormalVectors(2,iEdge) = yCell(cell2) - yCell(cell1)
edgeNormalVectors(3,iEdge) = zCell(cell2) - zCell(cell1)
- call unit_vector_in_R3(edgeNormalVectors(:,iEdge))
+ call unit_vec_in_R3(edgeNormalVectors(:,iEdge))
edgeLocations(1,iEdge) = 0.5*(xCell(cell2) + xCell(cell1))
edgeLocations(2,iEdge) = 0.5*(yCell(cell2) + yCell(cell1))
edgeLocations(3,iEdge) = 0.5*(zCell(cell2) + zCell(cell1))
@@ -166,10 +169,10 @@
rHat = radialUnitVectors(:,iCell)
normalDotRHat = sum(edgeNormalVectors(:,iEdge)*rHat)
xHatPlane = edgeNormalVectors(:,iEdge) - normalDotRHat*rHat
- call unit_vector_in_R3(xHatPlane)
+ call unit_vec_in_R3(xHatPlane)
call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
- call unit_vector_in_R3(yHatPlane) ! just to be sure...
+ call unit_vec_in_R3(yHatPlane) ! just to be sure...
cellTangentPlane(:,1,iCell) = xHatPlane
cellTangentPlane(:,2,iCell) = yHatPlane
end do
@@ -194,7 +197,7 @@
! coefficients - the coefficients needed to perform interpolation of the funciton
! at destination points yet to be specified
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_location_2D_scalar_constant_computeCoefficients(pointCount, coeffCount, &
+ subroutine rbfInterp_loc_2D_sca_const_compCoeffs(pointCount, coeffCount, &
points, fieldValues, alpha, coefficients)
integer, intent(in) :: pointCount, coeffCount
@@ -231,7 +234,7 @@
call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
coefficients(1:matrixSize), pivotIndices(1:matrixSize))
- end subroutine rbfInterp_location_2D_scalar_constant_computeCoefficients
+ end subroutine rbfInterp_loc_2D_sca_const_compCoeffs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: Compute interpolation coefficients in 2D that can be used to
@@ -252,7 +255,7 @@
! coefficients - the coefficients needed to perform interpolation of the funciton
! at destination points yet to be specified
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_location_2D_scalar_linear_computeCoefficients(pointCount, coeffCount, &
+ subroutine rbfInterp_loc_2D_sca_lin_compCoeffs(pointCount, coeffCount, &
points, fieldValues, alpha, coefficients)
integer, intent(in) :: pointCount, coeffCount
@@ -290,11 +293,11 @@
call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
coefficients(1:matrixSize), pivotIndices(1:matrixSize))
- end subroutine rbfInterp_location_2D_scalar_linear_computeCoefficients
+ end subroutine rbfInterp_loc_2D_sca_lin_compCoeffs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: Evalute a scalar function in 2D using coefficients computed in
- ! rbfInterp_location_2D_scalar_constant_computeCoefficients. This
+ ! rbfInterp_loc_2D_sca_const_compCoeffs. This
! function can be called repeatedly with different destination points
! to quickly evaluate the interpolating function using the same
! coefficients. This is useful for finding the location on the the
@@ -318,13 +321,12 @@
! derivs - the value of the function, the 2 components of its Jacobian and
! the 3 unique components of its Hessian at the evaluationPoint
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_location_2D_scalar_constant_evaluateWithDerivs(fieldCount, coeffCount, &
+ subroutine rbfInterp_loc_2D_sca_const_evalWithDerivs(fieldCount, coeffCount, &
pointCount, coefficients, evaluationPoint, points, alpha, derivs)
integer, intent(in) :: fieldCount, coeffCount, pointCount
real(kind=RKIND), dimension(coeffCount, fieldCount), intent(in) :: coefficients
real(kind=RKIND), dimension(2), intent(in) :: evaluationPoint
real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
- integer, intent(in) :: polynomialFlag
real(kind=RKIND), intent(in) :: alpha
real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
@@ -359,11 +361,11 @@
end if
end do
derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:)
- end subroutine rbfInterp_location_2D_scalar_constant_evaluateWithDerivs
+ end subroutine rbfInterp_loc_2D_sca_const_evalWithDerivs
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: Evalute a scalar function in 2D using coefficients computed in
- ! rbfInterp_location_2D_scalar_constant_computeCoefficients. This
+ ! rbfInterp_loc_2D_sca_const_compCoeffs. This
! function can be called repeatedly with different destination points
! to quickly evaluate the interpolating function using the same
! coefficients. This is useful for finding the location on the the
@@ -387,13 +389,12 @@
! derivs - the value of the function, the 2 components of its Jacobian and
! the 3 unique components of its Hessian at the evaluationPoint
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_location_2D_scalar_linear_evaluateWithDerivs(fieldCount, coeffCount, &
+ subroutine rbfInterp_loc_2D_sca_lin_evalWithDerivs(fieldCount, coeffCount, &
pointCount, coefficients, evaluationPoint, points, alpha, derivs)
integer, intent(in) :: fieldCount, coeffCount, pointCount
real(kind=RKIND), dimension(coeffCount, fieldCount), intent(in) :: coefficients
real(kind=RKIND), dimension(2), intent(in) :: evaluationPoint
real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
- integer, intent(in) :: polynomialFlag
real(kind=RKIND), intent(in) :: alpha
real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
@@ -433,10 +434,270 @@
derivs(2,:) = derivs(2,:) + coefficients(pointCount+2,:)
derivs(3,:) = derivs(3,:) + coefficients(pointCount+3,:)
- end subroutine rbfInterp_location_2D_scalar_linear_evaluateWithDerivs
+ end subroutine rbfInterp_loc_2D_sca_lin_evalWithDerivs
- subroutine rbfInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
- sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of scalar 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 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 boundary
+ ! conditions (or no boundaries). 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
+ ! 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_sca_const_dir_compCoeffs( &
+ pointCount, sourcePoints, destinationPoint, alpha, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount), intent(out) :: coefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+1 !! 1 extra space for constant
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ rhs = 0.0
+ coeffs = 0.0
+
+ call setUpScalarRBFDirichletMatrixAndRHS(pointCount, sourcePoints, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0
+ dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
+ end do
+
+ rhs(pointCount+1) = 1.0
+
+ ! solve each linear system
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine rbfInterp_func_3D_sca_const_dir_compCoeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in a plane in 3D that can be used to
+ ! interpolate a number of scalar 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 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 both Dirichlet (or no)
+ ! boundary conditions. The interpolation is performed with basis functions that are
+ ! RBFs plus constant and linear. All points are projected into the plane given by the
+ ! planeBasisVectors.
+ ! 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 points will be projected into the plane given by
+ ! planeBasisVectors
+ ! destinationPoint - the point in 3D where the interpolation will be performed. The
+ ! destinationPoint will be 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_sca_lin_dir_compCoeffs( &
+ pointCount, sourcePoints, destinationPoint, &
+ alpha, planeBasisVectors, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ 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), intent(out) :: coefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ rhs = 0.0
+ coeffs = 0.0
+
+ call setUpScalarRBFDirichletMatrixAndRHS(pointCount, sourcePoints, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0
+ dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
+ dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
+ dirichletMatrix(pointCount+1:pointCount+3,i) &
+ = dirichletMatrix(i,pointCount+1:pointCount+3)
+ end do
+
+ rhs(pointCount+1) = 1.0
+ rhs(pointCount+2) = sum(destinationPoint(1:3)*planeBasisVectors(1,:))
+ rhs(pointCount+3) = sum(destinationPoint(1:3)*planeBasisVectors(2,:))
+
+ ! solve each linear system
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine rbfInterp_func_3DPlane_sca_lin_dir_compCoeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of scalar 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 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 both Dirichlet (or no)
+ ! boundary conditions. The interpolation is performed with basis functions that are
+ ! RBFs plus constant and linear.
+ ! 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
+ ! 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_sca_lin_dir_compCoeffs(pointCount, &
+ sourcePoints, destinationPoint, alpha, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount), intent(out) :: coefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ rhs = 0.0
+ coeffs = 0.0
+
+ call setUpScalarRBFDirichletMatrixAndRHS(pointCount, sourcePoints, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0
+ dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
+ dirichletMatrix(pointCount+1:pointCount+4,i) &
+ = dirichletMatrix(i,pointCount+1:pointCount+4)
+ end do
+
+ rhs(pointCount+1) = 1.0
+ rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+
+ ! solve each linear system
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine rbfInterp_func_3D_sca_lin_dir_compCoeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of scalar 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 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 both Dirichlet and Neumann
+ ! boundary conditions. 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
+ ! isInterface - a logical array indicating which of the source points (if any) are at
+ ! at the domain interface. These points and their normals will be used to compute the
+ ! neumannCoefficients below
+ ! interfaceNormals - a 3D normal vector for each sourcePoint. These vectors are only used
+ ! at points where isInterface == .true., and can take arbitrary values elsewehere. The
+ ! normal vector is used to compute coefficients for the normal derivative of the
+ ! interpolating function in order to impose the Neumann Boundary condition
+ ! 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:
+ ! dirichletCoefficients - the coefficients used to interpolate a function with Dirichlet
+ ! boundary conditions to the specified destinationPoint
+ ! neumannCoefficients - the coefficients used to interpolate a function with Neumann
+ ! boundary conditions to the specified destinationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rbfInterp_func_3D_sca_const_dirNeu_compCoeffs( &
+ pointCount, sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
integer, intent(in) :: pointCount
@@ -503,10 +764,47 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeScalarRBFAndConstantCoefficients
+ end subroutine rbfInterp_func_3D_sca_const_dirNeu_compCoeffs
- subroutine rbfInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
- sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in a plane in 3D that can be used to
+ ! interpolate a number of scalar 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 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 both Dirichlet and Neumann
+ ! boundary conditions. The interpolation is performed with basis functions that are
+ ! RBFs plus constant and linear. All points are projected into the plane given by the
+ ! planeBasisVectors.
+ ! 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 will be projected into the plane given by
+ ! planeBasisVectors
+ ! isInterface - a logical array indicating which of the source points (if any) are at
+ ! at the domain interface. These points and their normals will be used to compute the
+ ! neumannCoefficients below
+ ! interfaceNormals - a 3D normal vector for each sourcePoint. These vectors are only used
+ ! at points where isInterface == .true., and can take arbitrary values elsewehere. The
+ ! normal vector is used to compute coefficients for the normal derivative of the
+ ! interpolating function in order to impose the Neumann Boundary condition
+ ! destinationPoint - the point in 3D where the interpolation will be performed. The
+ ! destinationPoint will be 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:
+ ! dirichletCoefficients - the coefficients used to interpolate a function with Dirichlet
+ ! boundary conditions to the specified destinationPoint
+ ! neumannCoefficients - the coefficients used to interpolate a function with Neumann
+ ! boundary conditions to the specified destinationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs( &
+ pointCount, sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients)
integer, intent(in) :: pointCount
@@ -583,9 +881,41 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeScalarRBFAndPlaneCoefficients
+ end subroutine rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs
- subroutine rbfInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 3D that can be used to
+ ! interpolate a number of scalar 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 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 both Dirichlet and Neumann
+ ! boundary conditions. The interpolation is performed with basis functions that are
+ ! RBFs plus constant and linear.
+ ! 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
+ ! isInterface - a logical array indicating which of the source points (if any) are at
+ ! at the domain interface. These points and their normals will be used to compute the
+ ! neumannCoefficients below
+ ! interfaceNormals - a 3D normal vector for each sourcePoint. These vectors are only used
+ ! at points where isInterface == .true., and can take arbitrary values elsewehere. The
+ ! normal vector is used to compute coefficients for the normal derivative of the
+ ! interpolating function in order to impose the Neumann Boundary condition
+ ! 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:
+ ! dirichletCoefficients - the coefficients used to interpolate a function with Dirichlet
+ ! boundary conditions to the specified destinationPoint
+ ! neumannCoefficients - the coefficients used to interpolate a function with Neumann
+ ! boundary conditions to the specified destinationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
@@ -659,9 +989,36 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeScalarRBFAndLinearCoefficients
+ end subroutine rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs
- subroutine rbfInterp_computeVectorDirichletRBFCoefficients(pointCount, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 boundary
+ ! conditions (or no boundaries). 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
+ ! 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.
+ ! 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_dir_compCoeffs(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients)
@@ -717,9 +1074,41 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeVectorDirichletRBFCoefficients
+ end subroutine rbfInterp_func_3D_vec_const_dir_compCoeffs
- subroutine rbfInterp_computeVectorDirichletPlanarRBFCoefficients(pointCount, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 boundary
+ ! conditions (or no boundaries). 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
+ ! 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. 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_dir_compCoeffs(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, planeBasisVectors, coefficients)
@@ -794,9 +1183,44 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeVectorDirichletPlanarRBFCoefficients
+ end subroutine rbfInterp_func_3DPlane_vec_const_dir_compCoeffs
- subroutine rbfInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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)
@@ -856,11 +1280,51 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeVectorFreeSlipRBFCoefficients
+ end subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs
- subroutine rbfInterp_computeVectorFreeSlipPlanarRBFCoefficients(pointCount, &
- sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
- alpha, planeBasisVectors, coefficients)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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
@@ -937,7 +1401,7 @@
deallocate(coeffs)
deallocate(pivotIndices)
- end subroutine rbfInterp_computeVectorFreeSlipPlanarRBFCoefficients
+ end subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
!!!!!!!!!!!!!!!!!!!!!
! private subroutines
@@ -980,6 +1444,35 @@
end subroutine evaluateRBFAndDerivs
+ subroutine setUpScalarRBFDirichletMatrixAndRHS(pointCount, sourcePoints, destinationPoint, &
+ alpha, dirichletMatrix, rhs)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: dirichletMatrix
+ real(kind=RKIND), dimension(pointCount), intent(out) :: rhs
+
+ integer :: i, j
+
+ real(kind=RKIND) :: rSquared, rbfValue
+
+ do j = 1, pointCount
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ dirichletMatrix(i,j) = rbfValue
+ end do
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j) = evaluateRBF(rSquared)
+ end do
+
+ end subroutine setUpScalarRBFDirichletMatrixAndRHS
+
subroutine setUpScalarRBFMatrixAndRHS(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletMatrix, neumannMatrix, rhs)
@@ -1105,13 +1598,13 @@
end subroutine setUpVectorFreeSlipRBFMatrixAndRHS
- subroutine unit_vector_in_R3(xin)
+ subroutine unit_vec_in_R3(xin)
implicit none
real (kind=RKIND), intent(inout) :: xin(3)
real (kind=RKIND) :: mag
mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
xin(:) = xin(:) / mag
- end subroutine unit_vector_in_R3
+ end subroutine unit_vec_in_R3
subroutine cross_product_in_R3(p_1,p_2,p_out)
real (kind=RKIND), intent(in) :: p_1 (3), p_2 (3)
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-05 03:51:01 UTC (rev 248)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-05-05 09:41:20 UTC (rev 249)
@@ -80,7 +80,6 @@
subroutine test_ibInterp_get2DRBF
integer, parameter :: pointCount = 16, coeffCount = 19, dimensions = 2, fieldCount = 1
- integer, parameter :: polynomialFlag = kMPAS_IBRBFPolynomialLinear
real(kind=RKIND), parameter :: alpha = 1.35, epsilon = 1e-6
real(kind=RKIND), dimension(pointCount,dimensions) :: points
real(kind=RKIND), dimension(pointCount) :: fieldValues
@@ -103,8 +102,8 @@
print *, "points:", points
print *, "fieldValues:", fieldValues
- call ibInterp_get2DRBFCoefficients(pointCount, coeffCount, dimensions, points, &
- fieldValues, polynomialFlag, alpha, coefficients)
+ call rbfInterp_loc_2D_sca_lin_compCoeffs(pointCount, coeffCount, points, &
+ fieldValues, alpha, coefficients)
print *, "coefficients:", coefficients
@@ -116,8 +115,8 @@
do xIndex = -1,1
evaluationPoint(pointIndex,1) = epsilon*xIndex + x0
evaluationPoint(pointIndex,2) = epsilon*yIndex + y0
- call ibInterp_evaluate2DRBFWithDerivs(fieldCount, coeffCount, pointCount, &
- coefficients, evaluationPoint(pointIndex,:), points, polynomialFlag, alpha, &
+ call rbfInterp_loc_2D_sca_lin_evalWithDerivs(fieldCount, coeffCount, pointCount, &
+ coefficients, evaluationPoint(pointIndex,:), points, alpha, &
derivs(pointIndex,:))
pointIndex = pointIndex + 1
end do
@@ -204,11 +203,11 @@
print *, "testType: ", testType
destinationPoint = sourcePoints(13,:)
if(testType == 1) then
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_const_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
else if(testType == 2) then
- call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
end if
@@ -221,11 +220,11 @@
destinationPoint = sourcePoints(32,:) + epsilon*interfaceNormals(1,:)
if(testType == 1) then
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_const_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
else if(testType == 2) then
- call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
end if
@@ -240,11 +239,11 @@
destinationPoint = sourcePoints(32,:) - epsilon*interfaceNormals(1,:)
if(testType == 1) then
- call ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_const_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
else if(testType == 2) then
- call ibInterp_computeScalarRBFAndLinearCoefficients(pointCount, &
+ call rbfInterp_func_3D_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, dirichletCoefficients, neumannCoefficients)
end if
@@ -314,7 +313,7 @@
planarBasisVectors(2,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
destinationPoint = sourcePoints(9,:)
- call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ call rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
@@ -325,7 +324,7 @@
destinationPoint = sourcePoints(12,:) + epsilon*interfaceNormals(1,:)
- call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ call rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
@@ -338,7 +337,7 @@
destinationPoint = sourcePoints(12,:) - epsilon*interfaceNormals(1,:)
- call ibInterp_computeScalarRBFAndPlaneCoefficients(pointCount, &
+ call rbfInterp_func_3DPlane_sca_lin_dirNeu_compCoeffs(pointCount, &
sourcePoints, isInterface, interfaceNormals, destinationPoint, &
alpha, planarBasisVectors, dirichletCoefficients, neumannCoefficients)
@@ -379,7 +378,7 @@
destinationPoint = sourcePoints(3,:)
- call ibInterp_computeVectorDirichletRBFCoefficients(pointCount, &
+ call rbfInterp_func_3D_vec_const_dir_compCoeffs(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients)
@@ -507,7 +506,7 @@
destinationPoint = sourcePoints(10,:)
- call ibInterp_computeVectorDirichletRBFCoefficients(pointCount, &
+ call rbfInterp_func_3D_vec_const_dir_compCoeffs(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients)
@@ -576,7 +575,7 @@
destinationPoint = sourcePoints(7,:)
- call ibInterp_computeVectorDirichletPlanarRBFCoefficients(pointCount, &
+ call rbfInterp_func_3DPlane_vec_const_dir_compCoeffs(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, planeBasisVectors, coefficients)
@@ -593,144 +592,5 @@
end subroutine test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
-! subroutine test_ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &
-! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
-! alpha, coefficients, dminfo)
-!
-! 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
-! 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 setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
-! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
-! alpha, matrix, rhs)
-!
-! 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
-! 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_computeVectorFreeSlipRBFCoefficients
-!
-! subroutine test_ibInterp_computeVectorFreeSlipPlanarRBFCoefficients(pointCount, &
-! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
-! alpha, planeBasisVectors, coefficients, dminfo)
-!
-! 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
-! type (dm_info), intent(in) :: dminfo
-!
-! 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
-! 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(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, rhs)
-!
-! 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
-! rhs(pointCount+1,:) = planeBasisVectors(1,:)
-! rhs(pointCount+2,:) = planeBasisVectors(2,:)
-!
-! ! solve each linear system
-! call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
-! call LEGS(matrix, 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(rhs)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine test_ibInterp_computeVectorFreeSlipPlanarRBFCoefficients
-!
-
end module interp_tests
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-05-05 03:51:01 UTC (rev 248)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-05-05 09:41:20 UTC (rev 249)
@@ -3,7 +3,7 @@
use grid_types
use configure
use constants
- use IB_interpolation
+ use RBF_interpolation
implicit none
@@ -277,7 +277,7 @@
tangentPlane(1,:) = cellTangentPlane(:,1,iCell)
tangentPlane(2,:) = cellTangentPlane(:,2,iCell)
- call ibInterp_computeVectorDirichletPlanarRBFCoefficients(pointCount, &
+ call rbfInterp_func_3DPlane_vec_const_dir_compCoeffs(pointCount, &
edgeOnCellLocations(1:pointCount,:), &
edgeOnCellNormals(1:pointCount,:), &
cellCenter, alpha, tangentPlane, coeffs(1:pointCount,:))
</font>
</pre>