<p><b>xylar@lanl.gov</b> 2010-04-21 16:12:14 -0600 (Wed, 21 Apr 2010)</p><p>Checking in the new interpolation module I'm working on just for safe keeping...<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F         (rev 0)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F        2010-04-21 22:12:14 UTC (rev 197)
@@ -0,0 +1,952 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+module IB_interpoation
+ use dmpar
+
+ implicit none
+ private
+ save
+
+ integer, parameter, public :: &
+ MPAS_i4 = selected_int_kind(6), &
+ MPAS_r8 = selected_real_kind(13)
+
+ public :: ibInterp_initialize, ibInterp_get2DRBFCoefficients, &
+ ibInterp_evaluate2DRBFWithDerivs, &
+ ibInterp_computeScalarRBFCoefficients, &
+ ibInterp_computeVectorNoSlipCoefficientsRBF, &
+ ibInterp_computeVectorFreeSlipCoefficientsRBF
+
+ integer(MPAS_i4), public, parameter :: &
+ kMPAS_IBRBFPolynomialConstant = 0, &
+ kMPAS_IBRBFPolynomialLinear = 1
+
+
+ contains
+
+ subroutine ibInterp_initialize(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose:
+ !
+ ! Input:
+ !
+ ! Output:
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (grid_meta), intent(inout) :: grid
+
+ integer(MPAS_i4) :: nCells, nEdges
+ integer(MPAS_i4), dimension(:,:), pointer :: cellsOnEdge
+ integer(MPAS_i4) :: iEdge, iCell, cell1, cell2
+ real(MPAS_r8), dimension(:), pointer :: xCell, yCell, zCell
+ real(MPAS_r8), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &
+ edgeLocations
+ real(MPAS_r8), dimension(3) :: xHatPlane, yHatPlane
+ real(MPAS_r8) :: normalDotRHat
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+
+ radialUnitVectors => grid % radialUnitVectors % array
+ edgeNormalVectors => grid % edgeNormalVectors % array
+ edgeLocations => grid % edgeLocations % array
+ cellTangentPlane => grid % cellTangentPlane % array
+
+ ! init arrays
+ edgeNormalVectors = 0
+ radialUnitVectors = 0
+
+ ! loop over all cells to be solved on this block
+ do iCell=1,nCells
+ radialUnitVectors(iCell,1) = xCell(iCell)
+ radialUnitVectors(iCell,2) = yCell(iCell)
+ radialUnitVectors(iCell,3) = zCell(iCell)
+ call unit_vector_in_R3(radialUnitVectors(iCell,:))
+ end do
+
+ do iEdge = 1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ edgeNormalVectors(iEdge,1) = xCell(cell2) - xCell(cell1)
+ edgeNormalVectors(iEdge,2) = yCell(cell2) - yCell(cell1)
+ edgeNormalVectors(iEdge,3) = zCell(cell2) - zCell(cell1)
+ edgeLocations(iEdge,1) = 0.5*(xCell(cell2) + xCell(cell1))
+ edgeLocations(iEdge,2) = 0.5*(yCell(cell2) + yCell(cell1))
+ edgeLocations(iEdge,3) = 0.5*(zCell(cell2) + zCell(cell1))
+ call unit_vector_in_R3(edgeNormalVectors(iEdge,:))
+ end do
+
+ do iCell=1,nCells
+ iEdge = edgesOnCell(1,iCell)
+ ! xHat and yHat are a local basis in the plane of the horizontal cell
+ ! we arbitrarily choose xHat to point toward the first edge
+ normalDotRHat = sum(edgeNormalVectors(iEdge,:)*radialUnitVectors(iCell,:))
+ xHatPlane = edgeNormalVectors(iEdge,:) - normalDotRHat*radialUnitVectors(iCell,1)
+ call unit_vector_in_R3(xHatPlane)
+
+ call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
+ call unit_vector_in_R3(yHatPlane) ! just to be sure...
+ cellTangentPlane(iCell, 1, :) = xHatPlane
+ cellTangentPlane(iCell, 2, :) = yHatPlane
+ end do
+
+ end subroutine ibInterp_initialize
+
+ subroutine ibInterp_getRBFCoefficients(pointCount, coeffCount, dimensions, points, fieldValues, &
+ polynomialFlag, alpha, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, coeffCount, dimensions
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: points
+ real(MPAS_r8), dimension(pointCount), intent(in) :: fieldValues
+ integer(MPAS_i4), intent(in) :: polynomialFlag
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(coeffCount), intent(out) :: coefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: i, j, matrixSize
+ real(MPAS_r8), dimension(pointCount+3,pointCount+3) :: matrix
+ real(MPAS_r8), dimension(pointCount+3) :: rhs
+ integer(MPAS_i4), dimension(pointCount+3) :: pivotIndices
+ real(MPAS_r8) :: rSquared
+
+ coefficients = 0.0_MPAS_r8
+ matrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+
+ rhs(1:pointCount) = fieldValues
+
+ do j=1,pointCount
+ do i=j,pointCount
+ rSquared = sum((points(j,:) - points(i,:))**2)/alpha**2
+ matrix(i,j) = evaluateRBF(rSquared)
+ matrix(j,i) = matrix(i,j)
+ end do
+ end do
+ if(polynomialFlag == kMPAS_IBRBFPolynomialConstant) then
+ matrixSize = pointCount+1
+ do j=1,pointCount
+ matrix(pointCount+1,j) = 1.0_MPAS_r8
+ matrix(j,pointCount+1) = 1.0_MPAS_r8
+ end do
+ else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
+ matrixSize = pointCount+3
+ do j=1,pointCount
+ matrix(pointCount+1,j) = 1.0_MPAS_r8
+ matrix(pointCount+2,j) = points(j,1)
+ matrix(pointCount+3,j) = points(j,2)
+ matrix(j,pointCount+1:pointCount+3) = matrix(pointCount+1:pointCount+3, j)
+ end do
+ else
+ call exit_MPAS("ibInterp_get2DRBFCoefficients: not a valid polynomialFlag", dminfo)
+ end if
+
+ call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
+ coefficients(1:matrixSize), pivotIndices(1:matrixSize))
+
+ end subroutine ibInterp_get2DRBFCoefficients
+
+ subroutine ibInterp_evaluate2DRBFWithDerivs(fieldCount, coeffCount, pointCount, &
+ coefficients, evaluationPoint, points, polynomialFlag, alpha, derivs, dminfo)
+ integer(MPAS_i4), intent(in) :: fieldCount, coeffCount, pointCount
+ real(MPAS_r8), dimension(coeffCount, fieldCount), intent(in) :: coefficients
+ real(MPAS_r8), dimension(2), intent(in) :: evaluationPoint
+ real(MPAS_r8), dimension(pointCount,2), intent(in) :: points
+ integer(MPAS_i4), intent(in) :: polynomialFlag
+ real(MPAS_r8), intent(in) :: alpha
+ type (dm_info), intent(in) :: dminfo
+
+ real(MPAS_r8), dimension(6,fieldCount), intent(out) :: derivs
+
+ integer(MPAS_i4) :: pointIndex
+ real(MPAS_r8) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+ derivs = 0.0_MPAS_r8
+ do pointIndex = 1, pointCount
+ x = (evaluationPoint(1) - points(pointIndex,1))
+ y = (evaluationPoint(2) - points(pointIndex,2))
+ rSquared = x**2 + y**2
+ call evaluateRBFAndDerivs(rSquared/alpha**2, rbfValue, rbfDerivOverR, rbfSecondDeriv)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ rbfSecondDeriv = rbfSecondDeriv/alpha**2
+ if(rSquared/alpha**2 < 1e-7) then
+ ! terms 2,3 and 5 are zero by radial symmetry of the RBF
+ derivs(1,:) = derivs(1,:) + coefficients(pointIndex,:)*rbfValue
+ derivs(4,:) = derivs(4,:) + coefficients(pointIndex,:)*rbfSecondDeriv
+ derivs(6,:) = derivs(6,:) + coefficients(pointIndex,:)*rbfSecondDeriv
+ else
+ call evaluateRBFAndDerivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
+ derivs(1,:) = derivs(1,:) + coefficients(pointIndex,:)*rbfValue
+ derivs(2,:) = derivs(2,:) + coefficients(pointIndex,:)*rbfDerivOverR*x
+ derivs(3,:) = derivs(3,:) + coefficients(pointIndex,:)*rbfDerivOverR*y
+ derivs(4,:) = derivs(4,:) + coefficients(pointIndex,:) &
+ * (rbfSecondDeriv*x**2 + rbfDerivOverR*y**2)/rSquared
+ derivs(5,:) = derivs(5,:) + coefficients(pointIndex,:) &
+ * (rbfSecondDeriv - rbfDerivOverR)*x*y/rSquared
+ derivs(6,:) = derivs(6,:) + coefficients(pointIndex,:) &
+ * (rbfSecondDeriv*y**2 + rbfDerivOverR*x**2)/rSquared
+ end if
+ end do
+ if(polynomialFlag == kMPAS_IBRBFPolynomialConstant) then
+ derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:)
+ else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
+ derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:) &
+ + coefficients(pointCount+2,:)*evaluationPoint(1) &
+ + coefficients(pointCount+3,:)*evaluationPoint(2)
+ derivs(2,:) = derivs(2,:) + coefficients(pointCount+2,:)
+ derivs(3,:) = derivs(3,:) + coefficients(pointCount+3,:)
+ else
+ call exit_MPAS("ibInterp_get2DRBFCoefficients: not a valid polynomialFlag", dminfo)
+ end if
+ end subroutine ibInterp_evaluate2DRBFWithDerivs
+
+ subroutine ibInterp_computeScalarRBFAndConstantCoefficients(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ polynomialFlag, alpha, dirichletCoefficients, neumannCoefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4), parameter :: matrixSize = pointCount+1 !! 1 extra space for constant
+
+
+ real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+ rhsCopy = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ call setUpScalarRBFMatrixAndRHS(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), &
+ neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+ else
+ neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
+ end if
+ dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
+ neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
+ end do
+
+ rhs(pointCount+1) = 1.0_MPAS_r8
+
+ ! solve each linear system
+ rhsCopy = rhs
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ end subroutine ibInterp_computeScalarRBFAndConstantCoefficients
+
+ subroutine ibInterp_computeScalarRBFAndPlane(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(2,3) :: planeBasisVectors
+ real(MPAS_r8), dimension(pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4), parameter :: matrixSize = pointCount+3 ! extra space for constant and plane
+
+ real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+ rhsCopy = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ call setUpScalarRBFMatrixAndRHS(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletMatrix(1:pointCount,1:pointCount), &
+ neumannMatrix(1:pointCount,1:pointCount), rhs(1:pointCount))
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+ dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVector(1,:))
+ dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVector(2,:))
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+ neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
+ else
+ neumannMatrix(i,pointCount+1:pointCount+4) &
+ = dirichletMatrix(i,pointCount+1:pointCount+4)
+ end if
+ dirichletMatrix(pointCount+1:pointCount+4,i) &
+ = dirichletMatrix(i,pointCount+1:pointCount+4)
+ neumannMatrix(pointCount+1:pointCount+4,i) &
+ = neumannMatrix(i,pointCount+1:pointCount+4)
+ end do
+
+ rhs(pointCount+1) = 1.0_MPAS_r8
+ rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+
+ ! solve each linear system
+ rhsCopy = rhs
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ end subroutine ibInterp_computeScalarRBFCoefficients
+
+ subroutine ibInterp_computeScalarRBFAndLinear(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(3,2) :: planeBasisVectors
+ real(MPAS_r8), dimension(pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: i, j
+ integer(MPAS_i4), parameter :: matrixSize = pointCount+4 ! extra space for constant and plane
+
+
+ real(MPAS_r8), dimension(matrixSize,matrixSize), pointer :: dirichletMatrix, neumannMatrix
+ real(MPAS_r8), dimension(matrixSize), pointer :: rhs, rhsCopy, coeffs
+ integer(MPAS_i4), dimension(matrixSize), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+ dirichletMatrix = 0.0_MPAS_r8
+ neumannMatrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+ rhsCopy = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ do j = 1, pointCount
+ if(isInterface(j)) cycle
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfValue
+ end do
+ end do
+
+ do j = 1, pointCount
+ if(.not. isInterface(j)) cycle
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalDotX = sum(interfaceNormals(j,:) &
+ * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+ end do
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j) = evaluateRBF(rSquared)
+ end do
+
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+ else
+ neumannMatrix(i,pointCount+1) = dirichletMatrix(i,pointCount+1)
+ end if
+ dirichletMatrix(pointCount+1,i) = dirichletMatrix(i,pointCount+1)
+ neumannMatrix(pointCount+1,i) = neumannMatrix(i,pointCount+1)
+ end do
+
+ rhs(pointCount+1) = 1.0_MPAS_r8
+ else if(polynomialFlag == kMPAS_IBRBFPolynomialLinear) then
+ do i = 1, pointCount
+ dirichletMatrix(i,pointCount+1) = 1.0_MPAS_r8
+ dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0_MPAS_r8
+ neumannMatrix(i,pointCount+2:pointCount+4) = interfaceNormals(i,1:3)
+ else
+ neumannMatrix(i,pointCount+1:pointCount+4) &
+ = dirichletMatrix(i,pointCount+1:pointCount+4)
+ end if
+ dirichletMatrix(pointCount+1:pointCount+4,i) &
+ = dirichletMatrix(i,pointCount+1:pointCount+4)
+ neumannMatrix(pointCount+1:pointCount+4,i) &
+ = neumannMatrix(i,pointCount+1:pointCount+4)
+ end do
+
+ rhs(pointCount+1) = 1.0_MPAS_r8
+ rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+ else
+ call exit_MPAS("ibInterp_computeScalarRBFCoefficients: not a valid polynomialFlag", dminfo)
+ end if
+
+
+ ! solve each linear system
+ rhsCopy = rhs
+ call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ end subroutine ibInterp_computeScalarRBFCoefficients
+
+ subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &
+ dimensions, sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: matrixSize, i, j
+
+ real(MPAS_r8), dimension(:,:), pointer :: matrix
+ real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+ matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,dimensions))
+ allocate(coeffs(matrixSize,dimensions))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ do j = 1, pointCount
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ rbfValue = evaluateRBF(rSquared)
+ matrix(j,i) = rbfValue*unitVectorDotProduct
+ end do
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
+
+ do i = 1, pointCount
+ matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+ matrix(pointCount+1:pointCount+dimensions,i) &
+ = matrix(i,pointCount+1:pointCount+dimensions)
+ end do
+ do i = 1, dimensions
+ rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ end do
+
+ ! solve each linear system
+ do i = 1, dimensions
+ call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ end do
+ coefficients = coeffs(1:pointCount,:)
+
+ end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF
+
+ subroutine ibInterp_computeVectorNoSlipCoefficientsRBF(pointCount, &
+ dimensions, sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: matrixSize, i, j
+
+ real(MPAS_r8), dimension(:,:), pointer :: matrix
+ real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+ matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,dimensions))
+ allocate(coeffs(matrixSize,dimensions))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ do j = 1, pointCount
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ rbfValue = evaluateRBF(rSquared)
+ matrix(j,i) = rbfValue*unitVectorDotProduct
+ end do
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
+
+ do i = 1, pointCount
+ matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+ matrix(pointCount+1:pointCount+dimensions,i) &
+ = matrix(i,pointCount+1:pointCount+dimensions)
+ end do
+ do i = 1, dimensions
+ rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ end do
+
+ ! solve each linear system
+ do i = 1, dimensions
+ call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ end do
+ coefficients = coeffs(1:pointCount,:)
+
+ end subroutine ibInterp_computeVectorNoSlipCoefficientsRBF
+
+ !!!!!! Not implemented yet, identical to No Slip !!!!!!!!!
+ subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF(pointCount, &
+ dimensions, sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients, dminfo)
+
+ integer(MPAS_i4), intent(in) :: pointCount, dimensions, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ real(MPAS_r8), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(MPAS_r8), dimension(dimensions), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount, dimensions), intent(out) :: coefficients
+ type (dm_info), intent(in) :: dminfo
+
+ integer(MPAS_i4) :: matrixSize, i, j
+
+ real(MPAS_r8), dimension(:,:), pointer :: matrix
+ real(MPAS_r8), dimension(:,:), pointer :: rhs, coeffs
+ integer(MPAS_i4), dimension(:), pointer :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, unitVectorDotProduct
+
+ matrixSize = pointCount + dimensions !! extra spaces for constant 2D flow
+
+ allocate(matrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize,dimensions))
+ allocate(coeffs(matrixSize,dimensions))
+ allocate(pivotIndices(matrixSize))
+
+ matrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+ coeffs = 0.0_MPAS_r8
+
+ do j = 1, pointCount
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ rbfValue = evaluateRBF(rSquared)
+ matrix(j,i) = rbfValue*unitVectorDotProduct
+ end do
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
+
+ do i = 1, pointCount
+ matrix(i,pointCount+1:pointCount+dimensions) = unitVector(i,:)
+ matrix(pointCount+1:pointCount+dimensions,i) &
+ = matrix(i,pointCount+1:pointCount+dimensions)
+ end do
+ do i = 1, dimensions
+ rhs(pointCount+pointCount+i,i) = 1.0_MPAS_r8 ! the unit vector in the ith direction
+ end do
+
+ ! solve each linear system
+ do i = 1, dimensions
+ call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ end do
+ coefficients = coeffs(1:pointCount,:)
+
+ end subroutine ibInterp_computeVectorFreeSlipCoefficientsRBF
+
+
+
+!!!!!!!!!!!!!!!!!!!!!
+! private subroutines
+!!!!!!!!!!!!!!!!!!!!!
+
+ subroutine exit_MPAS(string, dminfo)
+ char(len=*), intent(in) :: string
+ type (dm_info), intent(in) :: dminfo
+
+ write(0,*) string
+ call dmpar_abort(dminfo)
+ end subroutine exit_MPAS
+
+ function evaluateRBF(rSquared) result(rbfValue)
+ real(MPAS_r8), intent(in) :: rSquared
+ real(MPAS_r8) :: rbfValue
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+
+ end function evaluateRBF
+
+ subroutine evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ real(MPAS_r8), intent(in) :: rSquared
+ real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+ rbfDerivOverR = -rbfValue**3
+
+ end subroutine evaluateRBFAndDeriv
+
+ subroutine evaluateRBFAndDerivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
+ real(MPAS_r8), intent(in) :: rSquared
+ real(MPAS_r8), intent(out) :: rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+ rbfDerivOverR = -rbfValue**3
+ rbfSecondDeriv = (2*rSquared-1)*rbfValue**5
+
+ end subroutine evaluateRBFAndDerivs
+
+ subroutine setUpScalarRBFMatrixAndRHS(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletMatrix, neumannMatrix, rhs)
+
+ integer(MPAS_i4), intent(in) :: pointCount, polynomialFlag
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(MPAS_r8), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(MPAS_r8), dimension(3), intent(in) :: destinationPoint
+ real(MPAS_r8), intent(in) :: alpha
+ real(MPAS_r8), dimension(pointCount,pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+ real(MPAS_r8), dimension(pointCount), intent(out) :: rhs
+
+ integer(MPAS_i4) :: i, j
+
+
+ integer(MPAS_i4), dimension(matrixSize) :: pivotIndices
+
+ real(MPAS_r8) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+ dirichletMatrix = 0.0_MPAS_r8
+ neumannMatrix = 0.0_MPAS_r8
+ rhs = 0.0_MPAS_r8
+
+ do j = 1, pointCount
+ if(isInterface(j)) cycle
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfValue
+ end do
+ end do
+
+ do j = 1, pointCount
+ if(.not. isInterface(j)) cycle
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalDotX = sum(interfaceNormals(j,:) &
+ * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+ end do
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j) = evaluateRBF(rSquared)
+ end do
+
+ end subroutine setUpScalarRBFMatrixAndRHS
+
+ subroutine unit_vector_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
+
+!======================================================================
+! BEGINNING OF CROSS_PRODUCT_IN_R3
+!======================================================================
+ subroutine cross_product_in_R3(p_1,p_2,p_out)
+
+!-----------------------------------------------------------------------
+! PURPOSE: compute p_1 cross p_2 and place in p_out
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! intent(in)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(in) :: &
+ p_1 (3), &
+ p_2 (3)
+
+!-----------------------------------------------------------------------
+! intent(out)
+!-----------------------------------------------------------------------
+ real (kind=RKIND), intent(out) :: &
+ p_out (3)
+
+ p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
+ p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
+ p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
+
+ end subroutine cross_product_in_R3
+!======================================================================
+! END OF CROSS_PRODUCT_IN_R3
+!======================================================================
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! !
+! Please Note: !
+! !
+! (1) This computer program is written by Tao Pang in conjunction with !
+! his book, "An Introduction to Computational Physics," published !
+! by Cambridge University Press in 1997. !
+! !
+! (2) No warranties, express or implied, are made for this program. !
+! !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!PROGRAM EX43
+!
+!
+! An example of solving linear equation set A(N,N)*X(N) = B(N)
+! with the partial-pivoting Gaussian elimination scheme. The
+! numerical values are for the Wheatstone bridge example discussed
+! in Section 4.1 in the book with all resistances being 100 ohms
+! and the voltage 200 volts.
+!
+! IMPLICIT NONE
+! INTEGER, PARAMETER :: N=3
+! INTEGER :: I,J
+! INTEGER, DIMENSION (N) :: INDX
+! REAL, DIMENSION (N) :: X,B
+! REAL, DIMENSION (N,N) :: A
+! DATA B /200.0,0.0,0.0/, &
+! ((A(I,J), J=1,N),I=1,N) /100.0,100.0,100.0,-100.0, &
+! 300.0,-100.0,-100.0,-100.0, 300.0/
+!
+! CALL LEGS (A,N,B,X,INDX)
+!
+! WRITE (6, "(F16.8)") (X(I), I=1,N)
+!END PROGRAM EX43
+
+
+SUBROUTINE LEGS (A,N,B,X,INDX)
+!
+! Subroutine to solve the equation A(N,N)*X(N) = B(N) with the
+! partial-pivoting Gaussian elimination scheme.
+! Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ integer(MPAS_i4), INTENT (IN) :: N
+ integer(MPAS_i4) :: I,J
+ integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
+ real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N) :: A
+ real(MPAS_r8), INTENT (INOUT), DIMENSION (N) :: B
+ real(MPAS_r8), INTENT (OUT), DIMENSION (N) :: X
+!
+ CALL ELGS (A,N,INDX)
+!
+ DO I = 1, N-1
+ DO J = I+1, N
+ B(INDX(J)) = B(INDX(J))-A(INDX(J),I)*B(INDX(I))
+ END DO
+ END DO
+!
+ X(N) = B(INDX(N))/A(INDX(N),N)
+ DO I = N-1, 1, -1
+ X(I) = B(INDX(I))
+ DO J = I+1, N
+ X(I) = X(I)-A(INDX(I),J)*X(J)
+ END DO
+ X(I) = X(I)/A(INDX(I),I)
+ END DO
+!
+END SUBROUTINE LEGS
+!
+
+
+
+! Updated 10/24/2001.
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! !
+! Please Note: !
+! !
+! (1) This computer program is written by Tao Pang in conjunction with !
+! his book, "An Introduction to Computational Physics," published !
+! by Cambridge University Press in 1997. !
+! !
+! (2) No warranties, express or implied, are made for this program. !
+! !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+SUBROUTINE MIGS (A,N,X,INDX)
+!
+! Subroutine to invert matrix A(N,N) with the inverse stored
+! in X(N,N) in the output. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ integer(MPAS_i4), INTENT (IN) :: N
+ integer(MPAS_i4) :: I,J,K
+ integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
+ real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N):: A
+ real(MPAS_r8), INTENT (OUT), DIMENSION (N,N):: X
+ real(MPAS_r8), DIMENSION (N,N) :: B
+!
+ DO I = 1, N
+ DO J = 1, N
+ B(I,J) = 0.0
+ END DO
+ END DO
+ DO I = 1, N
+ B(I,I) = 1.0
+ END DO
+!
+ CALL ELGS (A,N,INDX)
+!
+ DO I = 1, N-1
+ DO J = I+1, N
+ DO K = 1, N
+ B(INDX(J),K) = B(INDX(J),K)-A(INDX(J),I)*B(INDX(I),K)
+ END DO
+ END DO
+ END DO
+!
+ DO I = 1, N
+ X(N,I) = B(INDX(N),I)/A(INDX(N),N)
+ DO J = N-1, 1, -1
+ X(J,I) = B(INDX(J),I)
+ DO K = J+1, N
+ X(J,I) = X(J,I)-A(INDX(J),K)*X(K,I)
+ END DO
+ X(J,I) = X(J,I)/A(INDX(J),J)
+ END DO
+ END DO
+END SUBROUTINE MIGS
+
+
+SUBROUTINE ELGS (A,N,INDX)
+!
+! Subroutine to perform the partial-pivoting Gaussian elimination.
+! A(N,N) is the original matrix in the input and transformed matrix
+! plus the pivoting element ratios below the diagonal in the output.
+! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
+!
+ IMPLICIT NONE
+ integer(MPAS_i4), INTENT (IN) :: N
+ integer(MPAS_i4) :: I,J,K,ITMP
+ integer(MPAS_i4), INTENT (OUT), DIMENSION (N) :: INDX
+ real(MPAS_r8) :: C1,PI,PI1,PJ
+ real(MPAS_r8), INTENT (INOUT), DIMENSION (N,N) :: A
+ real(MPAS_r8), DIMENSION (N) :: C
+!
+! Initialize the index
+!
+ DO I = 1, N
+ INDX(I) = I
+ END DO
+!
+! Find the rescaling factors, one from each row
+!
+ DO I = 1, N
+ C1= 0.0
+ DO J = 1, N
+ !C1 = AMAX1(C1,ABS(A(I,J)))
+ C1 = MAX(C1,ABS(A(I,J)))
+ END DO
+ C(I) = C1
+ END DO
+!
+! Search the pivoting (largest) element from each column
+!
+ DO J = 1, N-1
+ PI1 = 0.0
+ DO I = J, N
+ PI = ABS(A(INDX(I),J))/C(INDX(I))
+ IF (PI.GT.PI1) THEN
+ PI1 = PI
+ K = I
+ ENDIF
+ END DO
+!
+! Interchange the rows via INDX(N) to record pivoting order
+!
+ ITMP = INDX(J)
+ INDX(J) = INDX(K)
+ INDX(K) = ITMP
+ DO I = J+1, N
+ PJ = A(INDX(I),J)/A(INDX(J),J)
+!
+! Record pivoting ratios below the diagonal
+!
+ A(INDX(I),J) = PJ
+!
+! Modify other elements accordingly
+!
+ DO K = J+1, N
+ A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
+ END DO
+ END DO
+ END DO
+!
+END SUBROUTINE ELGS
+
+end module IB_interpoation
+
</font>
</pre>