<p><b>dwj07@fsu.edu</b> 2011-10-17 11:18:33 -0600 (Mon, 17 Oct 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Renaming of src/operators is complete.<br>
<br>
Tested building of each core with no issues.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/source_renaming/src/core_hyd_atmos/module_advection.F
===================================================================
--- branches/source_renaming/src/core_hyd_atmos/module_advection.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_hyd_atmos/module_advection.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -538,14 +538,14 @@
ata = matmul(at,a)
! if (m == n) then
-! call migs(a,n,b,indx)
+! call MIGS(a,n,b,indx)
! else
- call migs(atha,n,atha_inv,indx)
+ call MIGS(atha,n,atha_inv,indx)
b = matmul(atha_inv,ath)
-! call migs(ata,n,ata_inv,indx)
+! call MIGS(ata,n,ata_inv,indx)
! b = matmul(ata_inv,at)
! end if
b_out(1:n,1:m) = b(1:n,1:m)
@@ -575,7 +575,7 @@
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
-SUBROUTINE MIGS (A,N,X,INDX)
+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.
@@ -617,7 +617,7 @@
X(J,I) = X(J,I)/A(INDX(J),J)
END DO
END DO
-END SUBROUTINE MIGS
+END SUBROUTine MIGS
SUBROUTINE ELGS (A,N,INDX)
Modified: branches/source_renaming/src/core_hyd_atmos/module_mpas_core.F
===================================================================
--- branches/source_renaming/src/core_hyd_atmos/module_mpas_core.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_hyd_atmos/module_mpas_core.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -108,7 +108,7 @@
use grid_types
use advection
use time_integration
- use RBF_interpolation
+ use rbf_interpolation
use vector_reconstruction
implicit none
@@ -122,9 +122,9 @@
call compute_state_diagnostics(block % state % time_levs(1) % state, mesh)
call compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call initialize_advection_rk(mesh)
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
+ call mpas_rbf_interp_initialize(mesh)
+ call mpas_init_reconstruct(mesh)
+ call mpas_reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
block % diag % uReconstructX % array, &
block % diag % uReconstructY % array, &
block % diag % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_hyd_atmos/module_time_integration.F
===================================================================
--- branches/source_renaming/src/core_hyd_atmos/module_time_integration.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_hyd_atmos/module_time_integration.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -321,7 +321,7 @@
!
block => domain % blocklist
do while (associated(block))
- call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
+ call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % diag % uReconstructX % array, &
block % diag % uReconstructY % array, &
block % diag % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/source_renaming/src/core_ocean/mpas_ocn_advection.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_ocean/mpas_ocn_advection.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -593,14 +593,14 @@
ata = matmul(at,a)
! if (m == n) then
-! call migs(a,n,b,indx)
+! call MIGS(a,n,b,indx)
! else
- call migs(atha,n,atha_inv,indx)
+ call MIGS(atha,n,atha_inv,indx)
b = matmul(atha_inv,ath)
-! call migs(ata,n,ata_inv,indx)
+! call MIGS(ata,n,ata_inv,indx)
! b = matmul(ata_inv,at)
! end if
b_out(1:n,1:m) = b(1:n,1:m)
@@ -630,7 +630,7 @@
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
-SUBROUTINE MIGS (A,N,X,INDX)
+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.
@@ -672,7 +672,7 @@
X(J,I) = X(J,I)/A(INDX(J),J)
END DO
END DO
-END SUBROUTINE MIGS
+END SUBROUTine MIGS
SUBROUTINE ELGS (A,N,INDX)
Modified: branches/source_renaming/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/source_renaming/src/core_ocean/mpas_ocn_mpas_core.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_ocean/mpas_ocn_mpas_core.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -199,7 +199,7 @@
subroutine mpas_init_block(block, mesh, dt)!{{{
use grid_types
- use RBF_interpolation
+ use rbf_interpolation
use vector_reconstruction
implicit none
@@ -214,9 +214,9 @@
call compute_mesh_scaling(mesh)
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
+ call mpas_rbf_interp_initialize(mesh)
+ call mpas_init_reconstruct(mesh)
+ call mpas_reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
block % state % time_levs(1) % state % uReconstructX % array, &
block % state % time_levs(1) % state % uReconstructY % array, &
block % state % time_levs(1) % state % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_rk4.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -347,7 +347,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
+ call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
block % state % time_levs(2) % state % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_split.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_ocean/mpas_ocn_time_integration_split.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -1137,7 +1137,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
+ call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
block % state % time_levs(2) % state % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F
===================================================================
--- branches/source_renaming/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -165,10 +165,10 @@
! subroutine call.
tracersIn(1:maxLevelCell(iCell))=tracers(iTracer,1:maxLevelCell(iCell),iCell)
- call CubicSplineCoefficients(posZMidZLevel, &
+ call mpas_cubic_spline_coefficients(posZMidZLevel, &
tracersIn, maxLevelCell(iCell), tracer2ndDer)
- call InterpolateCubicSpline( &
+ call mpas_interpolate_cubic_spline( &
posZMidZLevel, tracersIn, tracer2ndDer, maxLevelCell(iCell), &
posZTopZLevel, tracersOut, maxLevelCell(iCell)-1 )
Modified: branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_sw/mpas_sw_mpas_core.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -115,7 +115,7 @@
use grid_types
use time_integration
- use RBF_interpolation
+ use rbf_interpolation
use vector_reconstruction
implicit none
@@ -128,9 +128,9 @@
call sw_compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call compute_mesh_scaling(mesh)
- call rbfInterp_initialize(mesh)
- call init_reconstruct(mesh)
- call reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
+ call mpas_rbf_interp_initialize(mesh)
+ call mpas_init_reconstruct(mesh)
+ call mpas_reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
block % state % time_levs(1) % state % uReconstructX % array, &
block % state % time_levs(1) % state % uReconstructY % array, &
block % state % time_levs(1) % state % uReconstructZ % array, &
Modified: branches/source_renaming/src/core_sw/mpas_sw_time_integration.F
===================================================================
--- branches/source_renaming/src/core_sw/mpas_sw_time_integration.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/core_sw/mpas_sw_time_integration.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -229,7 +229,7 @@
call sw_compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
- call reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
+ call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
block % state % time_levs(2) % state % uReconstructZ % array, &
Modified: branches/source_renaming/src/operators/Makefile
===================================================================
--- branches/source_renaming/src/operators/Makefile 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/operators/Makefile 2011-10-17 17:18:33 UTC (rev 1089)
@@ -1,15 +1,15 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_vector_reconstruction.o module_spline_interpolation.o
+OBJS = mpas_rbf_interpolation.o mpas_vector_reconstruction.o mpas_spline_interpolation.o
all: operators
operators: $(OBJS)
ar -ru libops.a $(OBJS)
-module_vector_reconstruction.o: module_RBF_interpolation.o
-module_RBF_interpolation.o:
-module_spline_interpolation:
+mpas_vector_reconstruction.o: mpas_rbf_interpolation.o
+mpas_rbf_interpolation.o:
+mpas_spline_interpolation:
clean:
$(RM) *.o *.mod *.f90 libops.a
Deleted: branches/source_renaming/src/operators/module_RBF_interpolation.F
===================================================================
--- branches/source_renaming/src/operators/module_RBF_interpolation.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/operators/module_RBF_interpolation.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -1,1824 +0,0 @@
-module RBF_interpolation
- use dmpar
- use grid_types
-
- implicit none
- private
- save
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Purpose: perform interpolation of scalar and vector functions in 2D
-! and 3D using Radial Basis Functions (RBFs).
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- ! Initialize the geometry that will be useful from interpolation
- public :: rbfInterp_initialize
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Routines for perofrming interpolation in 2D (including Jacobian and Hessian)
- ! at locations that vary using a function that is fixed. This is useful
- ! 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_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
- ! in 3D) with support for either constant or constant and linear basis
- ! functions in addition to RBFs. In constrast to the two subroutines
- ! 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.
- ! 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(functionOrDerivAtSources*neumannCoefficients)
- ! where functionOrDerivAtSource = functionAtSources where .not.(isInterface)
- ! = functionNormalDerivAtSources where isInterface
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- 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
- ! in 3D) with support for only constant basis functions in addition to RBFs.
- ! (Linear basis functions can cause problems: a linear vortex flow u = y xHat - x yHat
- ! cannot be reconstructed from vectors that are only normal to the cell edges in 2D.
- ! Therefore, we don't support them). As with the scalar 3D subroutines
- ! 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.
- ! The user supplies to these routines a set of sourcePoints and unitVectors
- ! as well as a destinationPoint and, for the last 2 routines, flags for
- ! which points are tangent to the interface and which of the supplied unitVectors
- ! is the normal at the corresponding point.
- !
- ! The first two routines compute coefficients that can be used to interpolate
- ! a vector function to the destination point given "function dot unitVector" values
- ! at each source point. These routines are useful, for example, for reconstructing
- ! the full vector velocity at cell centers from the normal components of the velocity
- ! at cell faces (or cell edges in 2D), or for computing the full velocity at an
- ! immersed boundary image point based on the normal velocity at several faces and
- ! the full velocity at boundary points (e.g., a no-slip boundary condition).
- !
- ! The last two routines compute coefficients that can be used to interpolate
- ! a vector function to the destination point given "function dot unitVector" values
- ! at non-tangent source point and "dFunction/dn dot unitVector" values at
- ! tangent source point. These routines are useful, for example, for computing the
- ! full velocity at an immersed boundary image point based on the normal velocity at
- ! several faces, the normal velocity at boundary points (e.g., u dot n = 0 for a
- ! no-penetration boundary condition on a fixed boundary), and the normal derivative
- ! of the tangential components of velocity at the boundary points (e.g., a free-slip
- ! boundary condition).
- ! Pseudocode for function reconstruction at the destinationPoint is as follows
- ! dirichlet: functionAtDestination_i = sum_j(functionDotUnitVectorAtSources_j*coefficients_j,i)
- ! for i = x,y,z
- ! tangentNeumann: functionAtDestination_i &
- ! = sum_(j where .not isTangent) (functionDotUnitVectorAtSources_j*coefficients_j,i) &
- ! + sum_(j where isTangent) ((dFunction/dn dot UnitVector)_j*coefficients_j,i)
- ! for i = x,y,z
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- 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
-
- subroutine rbfInterp_initialize(grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: compute geometric fields that will be potentially useful for calling
- ! the interpolation routines
- !
- ! Input: the grid
- !
- ! Output:
- ! edgeNormalVectors - the unit vector at the center of each edge tangent to the sphere
- ! cellTangentPlane - 2 orthogonal unit vectors in the tangent plane of each cell
- ! The first unit vector is chosen to point toward the center of the first
- ! edge on the cell.
- ! localVerticalUnitVectors - the unit normal vector of the tangent plane at the center
- ! of each cell
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
-
- integer :: nCells, nEdges
- integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
- integer :: iEdge, iCell, cell1, cell2
- real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
- real(kind=RKIND), dimension(:,:), pointer :: localVerticalUnitVectors, edgeNormalVectors
- real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
- real(kind=RKIND), dimension(3) :: xHatPlane, yHatPlane, rHat
- real(kind=RKIND) :: normalDotRHat
- logical :: on_a_sphere
-
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- xEdge => grid % xEdge % array
- yEdge => grid % yEdge % array
- zEdge => grid % zEdge % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nCells = grid % nCells
- nEdges = grid % nEdges
- on_a_sphere = grid % on_a_sphere
-
- localVerticalUnitVectors => grid % localVerticalUnitVectors % array
- edgeNormalVectors => grid % edgeNormalVectors % array
- cellTangentPlane => grid % cellTangentPlane % array
-
- ! init arrays
- edgeNormalVectors = 0
- localVerticalUnitVectors = 0
-
- ! loop over all cells to be solved on this block
- do iCell=1,nCells
- if(on_a_sphere) then
- localVerticalUnitVectors(1,iCell) = xCell(iCell)
- localVerticalUnitVectors(2,iCell) = yCell(iCell)
- localVerticalUnitVectors(3,iCell) = zCell(iCell)
- call unit_vec_in_R3(localVerticalUnitVectors(:,iCell))
- else ! on a plane
- localVerticalUnitVectors(:,iCell) = (/ 0., 0., 1. /)
- end if
- end do
-
- do iEdge = 1,nEdges
- iCell = cellsOnEdge(1,iEdge) ! the normal vector points from the first cell toward the edge
- if(iCell == nCells+1) then ! this is a boundary edge
- ! the first cell bordering this edge is not real, use the second cell
- ! The normal should always point outward at boundaries, away from the valid cell center
- iCell = cellsOnEdge(2,iEdge)
- end if
- ! the normal points from the cell location to the edge location
- edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(iCell)
- edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(iCell)
- edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(iCell)
- call unit_vec_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
- rHat = localVerticalUnitVectors(:,iCell)
- normalDotRHat = sum(edgeNormalVectors(:,iEdge)*rHat)
- xHatPlane = edgeNormalVectors(:,iEdge) - normalDotRHat*rHat
- call unit_vec_in_R3(xHatPlane)
-
- call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
- call unit_vec_in_R3(yHatPlane) ! just to be sure...
- cellTangentPlane(:,1,iCell) = xHatPlane
- cellTangentPlane(:,2,iCell) = yHatPlane
- end do
-
- end subroutine rbfInterp_initialize
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: Compute interpolation coefficients in 2D that can be used to
- ! reconstruct a given scalar function at varying locations. This is useful
- ! 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.
- ! The reconstruction is performed with basis functions that are RBFs and constant
- ! Input:
- ! pointCount - the number of "source" points and functionValues supplied
- ! coeffCount - the size of coefficients, must be at least pointCount + 1
- ! points - the location of the "source" points in the 2D space where the values of
- ! the function are known
- ! fieldValues - the values of the function of interest at the points
- ! 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 needed to perform interpolation of the funciton
- ! at destination points yet to be specified
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_loc_2D_sca_const_compCoeffs(pointCount, coeffCount, &
- points, fieldValues, alpha, coefficients)
-
- integer, intent(in) :: pointCount, coeffCount
- real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
- real(kind=RKIND), dimension(pointCount), intent(in) :: fieldValues
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(coeffCount), intent(out) :: coefficients
-
- integer :: i, j, matrixSize
- real(kind=RKIND), dimension(pointCount+1,pointCount+1) :: matrix
- real(kind=RKIND), dimension(pointCount+1) :: rhs
- integer, dimension(pointCount+1) :: pivotIndices
- real(kind=RKIND) :: rSquared
-
- matrixSize = pointCount+1
- coefficients = 0.0
- matrix = 0.0
- rhs = 0.0
-
- 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
- do j=1,pointCount
- matrix(pointCount+1,j) = 1.0
- matrix(j,pointCount+1) = 1.0
- end do
-
- call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
- coefficients(1:matrixSize), pivotIndices(1:matrixSize))
-
- end subroutine rbfInterp_loc_2D_sca_const_compCoeffs
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: Compute interpolation coefficients in 2D that can be used to
- ! reconstruct a given scalar function at varying locations. This is useful
- ! 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.
- ! The reconstruction is performed with basis functions that are RBFs plus constant
- ! and linear
- ! Input:
- ! pointCount - the number of "source" points and functionValues supplied
- ! coeffCount - the size of coefficients, must be at least pointCount + 3
- ! points - the location of the "source" points in the 2D space where the values of
- ! the function are known
- ! fieldValues - the values of the function of interest at the points
- ! 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 needed to perform interpolation of the funciton
- ! at destination points yet to be specified
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_loc_2D_sca_lin_compCoeffs(pointCount, coeffCount, &
- points, fieldValues, alpha, coefficients)
-
- integer, intent(in) :: pointCount, coeffCount
- real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
- real(kind=RKIND), dimension(pointCount), intent(in) :: fieldValues
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(coeffCount), intent(out) :: coefficients
-
- integer :: i, j, matrixSize
- real(kind=RKIND), dimension(pointCount+3,pointCount+3) :: matrix
- real(kind=RKIND), dimension(pointCount+3) :: rhs
- integer, dimension(pointCount+3) :: pivotIndices
- real(kind=RKIND) :: rSquared
-
- coefficients = 0.0
- matrix = 0.0
- rhs = 0.0
-
- 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
- matrixSize = pointCount+3
- do j=1,pointCount
- matrix(pointCount+1,j) = 1.0
- 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
- call LEGS(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
- coefficients(1:matrixSize), pivotIndices(1:matrixSize))
-
- end subroutine rbfInterp_loc_2D_sca_lin_compCoeffs
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: Evalute a scalar function in 2D using coefficients computed in
- ! 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
- ! RBF reconstruction of a function (e.g., a height field) that minimizes
- ! the distance to a point in 3D space. The reconstruction is performed
- ! with basis functions that are RBFs and constant
- ! Input:
- ! fieldCount - the number fields to be evaluated. This is useful for reconstructing,
- ! for example, the x-, y- and z-components of a vector field at the same
- ! point in 2D
- ! coeffCount - the size of coefficients, must be at least pointCount + 1
- ! pointCount - the number of "source" points and functionValues supplied
- ! coefficients - the coefficients needed to perform interpolation of the funciton
- ! at the evaluationPoint
- ! evaluationPoint - the point in 2D where the function is to be reconstructed
- ! points - the location of the "source" points in the 2D space where the values of
- ! the function are known
- ! alpha - a constant that give the characteristic length scale of the RBFs,
- ! should be on the order of the distance between points
- ! Output:
- ! 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_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
- real(kind=RKIND), intent(in) :: alpha
-
- real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
-
- integer :: pointIndex
- real(kind=RKIND) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
-
- derivs = 0.0
- 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
- derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:)
- end subroutine rbfInterp_loc_2D_sca_const_evalWithDerivs
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: Evalute a scalar function in 2D using coefficients computed in
- ! 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
- ! RBF reconstruction of a function (e.g., a height field) that minimizes
- ! the distance to a point in 3D space. The reconstruction is performed
- ! with basis functions that are RBFs, constant and linear
- ! Input:
- ! fieldCount - the number fields to be evaluated. This is useful for reconstructing,
- ! for example, the x-, y- and z-components of a vector field at the same
- ! point in 2D
- ! coeffCount - the size of coefficients, must be at least pointCount + 1
- ! pointCount - the number of "source" points and functionValues supplied
- ! coefficients - the coefficients needed to perform interpolation of the funciton
- ! at the evaluationPoint
- ! evaluationPoint - the point in 2D where the function is to be reconstructed
- ! points - the location of the "source" points in the 2D space where the values of
- ! the function are known
- ! alpha - a constant that give the characteristic length scale of the RBFs,
- ! should be on the order of the distance between points
- ! Output:
- ! 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_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
- real(kind=RKIND), intent(in) :: alpha
-
- real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
-
- integer :: pointIndex
- real(kind=RKIND) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
-
- derivs = 0.0
- 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
- 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,:)
-
- end subroutine rbfInterp_loc_2D_sca_lin_evalWithDerivs
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! 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
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isInterface
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount), intent(out) :: &
- dirichletCoefficients, neumannCoefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
- real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
- integer, dimension(:), pointer :: pivotIndices
-
- matrixSize = pointCount+1 !! 1 extra space for constant
-
- allocate(dirichletMatrix(matrixSize,matrixSize))
- allocate(neumannMatrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize))
- allocate(rhsCopy(matrixSize))
- allocate(coeffs(matrixSize))
- allocate(pivotIndices(matrixSize))
-
- dirichletMatrix = 0.0
- neumannMatrix = 0.0
- rhs = 0.0
- rhsCopy = 0.0
- coeffs = 0.0
-
- 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
- if(isInterface(i)) then
- neumannMatrix(i,pointCount+1) = 0.0
- 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
-
- ! 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)
-
- deallocate(dirichletMatrix)
- deallocate(neumannMatrix)
- deallocate(rhs)
- deallocate(rhsCopy)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- end subroutine rbfInterp_func_3D_sca_const_dirNeu_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 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
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isInterface
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
- 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) :: &
- dirichletCoefficients, neumannCoefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
- real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
- integer, dimension(:), pointer :: pivotIndices
-
- matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
-
- allocate(dirichletMatrix(matrixSize,matrixSize))
- allocate(neumannMatrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize))
- allocate(rhsCopy(matrixSize))
- allocate(coeffs(matrixSize))
- allocate(pivotIndices(matrixSize))
-
- dirichletMatrix = 0.0
- neumannMatrix = 0.0
- rhs = 0.0
- rhsCopy = 0.0
- coeffs = 0.0
-
- 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
- dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
- dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
- if(isInterface(i)) then
- neumannMatrix(i,pointCount+1) = 0.0
- neumannMatrix(i,pointCount+2) = sum(interfaceNormals(i,1:3)*planeBasisVectors(1,:))
- neumannMatrix(i,pointCount+3) = sum(interfaceNormals(i,1:3)*planeBasisVectors(2,:))
- else
- neumannMatrix(i,pointCount+1:pointCount+3) &
- = dirichletMatrix(i,pointCount+1:pointCount+3)
- end if
- dirichletMatrix(pointCount+1:pointCount+3,i) &
- = dirichletMatrix(i,pointCount+1:pointCount+3)
- neumannMatrix(pointCount+1:pointCount+3,i) &
- = neumannMatrix(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
- rhsCopy = rhs
- call LEGS(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
- dirichletCoefficients = coeffs(1:pointCount)
-
- call LEGS(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
- neumannCoefficients = coeffs(1:pointCount)
-
- deallocate(dirichletMatrix)
- deallocate(neumannMatrix)
- deallocate(rhs)
- deallocate(rhsCopy)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- end subroutine rbfInterp_func_3DPlane_sca_lin_dirNeu_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 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)
-
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isInterface
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount), intent(out) :: &
- dirichletCoefficients, neumannCoefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
- real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
- integer, dimension(:), pointer :: pivotIndices
-
- matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
-
- allocate(dirichletMatrix(matrixSize,matrixSize))
- allocate(neumannMatrix(matrixSize,matrixSize))
- allocate(rhs(matrixSize))
- allocate(rhsCopy(matrixSize))
- allocate(coeffs(matrixSize))
- allocate(pivotIndices(matrixSize))
-
- dirichletMatrix = 0.0
- neumannMatrix = 0.0
- rhs = 0.0
- rhsCopy = 0.0
- coeffs = 0.0
-
- 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
- dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
- if(isInterface(i)) then
- neumannMatrix(i,pointCount+1) = 0.0
- 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
- 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)
-
- deallocate(dirichletMatrix)
- deallocate(neumannMatrix)
- deallocate(rhs)
- deallocate(rhsCopy)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- end subroutine rbfInterp_func_3D_sca_lin_dirNeu_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 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)
-
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-
- 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 setUpVectorDirichletRBFMatrixAndRHS(pointCount, 3, &
- sourcePoints, unitVectors, destinationPoint, &
- alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-
- do i = 1, pointCount
- matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
- matrix(pointCount+1:pointCount+3,i) &
- = matrix(i,pointCount+1:pointCount+3)
- end do
- do i = 1, 3
- rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
- end do
-
- ! solve each linear system
- do i = 1, 3
- 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_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 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)
-
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(2,3) :: planeBasisVectors
- real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
-
- 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 setUpVectorDirichletRBFMatrixAndRHS(pointCount, 2, &
- planarSourcePoints, planarUnitVectors, planarDestinationPoint, &
- alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-
- do i = 1, pointCount
- matrix(i,pointCount+1:pointCount+2) = planarUnitVectors(i,:)
- 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)
-
-
- do i = 1,3
- coefficients(:,i) = planeBasisVectors(1,i)*coeffs(1:pointCount,1) &
- + planeBasisVectors(2,i)*coeffs(1:pointCount,2)
- end do
-
- deallocate(matrix)
- deallocate(matrixCopy)
- deallocate(rhs)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- 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
- 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
-!!!!!!!!!!!!!!!!!!!!!
-
- function evaluateRBF(rSquared) result(rbfValue)
- real(kind=RKIND), intent(in) :: rSquared
- real(kind=RKIND) :: rbfValue
-
- ! inverse multiquadratic
- rbfValue = 1/sqrt(1 + rSquared)
-
- end function evaluateRBF
-
- subroutine evaluateRBFAndDeriv(rSquared, rbfValue, rbfDerivOverR)
- real(kind=RKIND), intent(in) :: rSquared
- real(kind=RKIND), intent(out) :: rbfValue, rbfDerivOverR
-
- ! inverse multiquadratic
- rbfValue = 1/sqrt(1 + rSquared)
- rbfDerivOverR = -rbfValue**3
-
- end subroutine evaluateRBFAndDeriv
-
- subroutine evaluateRBFAndDerivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
- real(kind=RKIND), intent(in) :: rSquared
- real(kind=RKIND), 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 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)
-
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isInterface
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: &
- dirichletMatrix, neumannMatrix
- real(kind=RKIND), dimension(pointCount), intent(out) :: rhs
-
- integer :: i, j
-
- real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
-
- do j = 1, pointCount
- if(isInterface(j)) then
- 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
- else
- 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 if
- end do
-
- do j = 1, pointCount
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j) = evaluateRBF(rSquared)
- end do
-
- end subroutine setUpScalarRBFMatrixAndRHS
-
- subroutine setUpVectorDirichletRBFMatrixAndRHS(pointCount, dimensions, &
- sourcePoints, unitVectors, destinationPoint, &
- alpha, matrix, rhs)
-
- integer, intent(in) :: pointCount, dimensions
- real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
- real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
- real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
- real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
-
- integer :: i, j
-
- real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
-
- do j = 1, pointCount
- do i = j, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- rbfValue = evaluateRBF(rSquared)
- unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- matrix(i,j) = rbfValue*unitVectorDotProduct
- matrix(j,i) = matrix(i,j)
- end do
- end do
-
- do j = 1, pointCount
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
- end do
-
- end subroutine setUpVectorDirichletRBFMatrixAndRHS
-
- subroutine setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, dimensions, &
- sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
- alpha, matrix, rhs)
-
- integer, intent(in) :: pointCount, dimensions
- real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isTangentToInterface
- integer, dimension(pointCount), intent(in) :: normalVectorIndex
- real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
- real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
- real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
-
- integer :: i, j
-
- real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalVector(dimensions), &
- normalDotX, unitVectorDotProduct
-
- do j = 1, pointCount
- if(isTangentToInterface(j)) then
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- 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
- else
- do i = 1, pointCount
- rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
- rbfValue = evaluateRBF(rSquared)
- unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
- matrix(i,j) = rbfValue*unitVectorDotProduct
- end do
- end if
- end do
-
- do j = 1, pointCount
- rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
- rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
- end do
-
- end subroutine setUpVectorFreeSlipRBFMatrixAndRHS
-
- 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_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)
- 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
-
-! 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, INTENT (IN) :: N
- integer :: I,J
- integer, INTENT (OUT), DIMENSION (N) :: INDX
- real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
- real(kind=RKIND), INTENT (INOUT), DIMENSION (N) :: B
- real(kind=RKIND), 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, INTENT (IN) :: N
- integer :: I,J,K
- integer, INTENT (OUT), DIMENSION (N) :: INDX
- real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
- real(kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
- real(kind=RKIND), 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, INTENT (IN) :: N
- integer :: I,J,K,ITMP
- integer, INTENT (OUT), DIMENSION (N) :: INDX
- real(kind=RKIND) :: C1,PI,PI1,PJ
- real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
- real(kind=RKIND), 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 RBF_interpolation
-
Deleted: branches/source_renaming/src/operators/module_spline_interpolation.F
===================================================================
--- branches/source_renaming/src/operators/module_spline_interpolation.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/operators/module_spline_interpolation.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -1,427 +0,0 @@
-module spline_interpolation
-
- implicit none
-
- private
-
- public :: CubicSplineCoefficients, InterpolateCubicSpline, &
- IntegrateCubicSpline, IntegrateColumnCubicSpline, InterpolateLinear, &
- TestInterpolate
-
-! Short Descriptions:
-
-! CubicSplineCoefficients: Compute second derivatives at nodes.
-! This must be run before any of the other cubic spine functions.
-
-! InterpolateCubicSpline: Compute cubic spline interpolation.
-
-! IntegrateCubicSpline: Compute a single integral from spline data.
-
-! IntegrateColumnCubicSpline: Compute multiple integrals from spline data.
-
-! InterpolateLinear: Compute linear interpolation.
-
-! TestInterpolate: Test spline interpolation subroutines.
-
- contains
-
- subroutine CubicSplineCoefficients(x,y,n,y2ndDer)
-
-! Given arrays x(1:n) and y(1:n) containing a function,
-! i.e., y(i) = f(x(i)), with x monotonically increasing
-! this routine returns an array y2ndDer(1:n) that contains
-! the second derivatives of the interpolating function at x(1:n).
-! This routine uses boundary conditions for a natural spline,
-! with zero second derivative on that boundary.
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n ! number of nodes
- real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! location of nodes
- y ! value at nodes
-
-! OUTPUT PARAMETERS:
-
- real(kind=RKIND), intent(out), dimension(n) :: &
- y2ndDer ! dy^2/dx^2 at each node
-
-! local variables:
-
- integer :: i
- real(kind=RKIND) :: &
- temp,xRatio,a(n)
-
- y2ndDer(1)=0.0
- y2ndDer(n)=0.0
- a(1)=0.0
-
- do i=2,n-1
- xRatio=(x(i)-x(i-1))/(x(i+1)-x(i-1))
- temp=1.0/(2.0+xRatio*y2ndDer(i-1))
- y2ndDer(i)=temp*(xRatio-1.0)
- a(i) = temp*(6.0*((y(i+1)-y(i))/(x(i+1)-x(i)) &
- -(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1)) &
- -xRatio*a(i-1))
- enddo
-
- do i=n-1,1,-1
- y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
- enddo
-
- end subroutine CubicSplineCoefficients
-
-
- subroutine InterpolateCubicSpline( &
- x,y,y2ndDer,n, &
- xOut,yOut,nOut)
-
-! Given the arrays x(1:n) and y(1:n), which tabulate a function,
-! and given the array y2ndDer(1:n), which is the output from
-! CubicSplineCoefficients above, this routine returns the
-! cubic-spline interpolated values of yOut(1:nOut) at xOut(1:nOut).
-! This subroutine assumes that both x and xOut are monotonically
-! increasing, and that all values of xOut are within the first and
-! last values of x.
-
-! INPUT PARAMETERS:
-
- real (kind=RKIND), dimension(n), intent(in) :: &
- x, &! node location, input grid
- y, &! interpolation variable, input grid
- y2ndDer ! 2nd derivative of y at nodes
-
- real (kind=RKIND), dimension(nOut), intent(in) :: &
- xOut ! node location, output grid
-
- integer, intent(in) :: &
- n, &! number of nodes, input grid
- nOut ! number of nodes, output grid
-
-! OUTPUT PARAMETERS:
-
- real (kind=RKIND), dimension(nOut), intent(out) :: &
- yOut ! interpolation variable, output grid
-
-! local variables:
-
- integer :: &
- kIn, kOut ! counters
-
- real (kind=RKIND) :: &
- a, b, h
-
- kOut = 1
-
- kInLoop: do kIn = 1,n-1
-
- h = x(kIn+1)-x(kIn)
-
- do while(xOut(kOut) < x(kIn+1))
-
- a = (x(kIn+1)-xOut(kOut))/h
- b = (xOut(kOut)-x (kIn) )/h
- yOut(kOut) = a*y(kIn) + b*y(kIn+1) &
- + ((a**3-a)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
- *(h**2)/6.0
-
- kOut = kOut + 1
-
- if (kOut>nOut) exit kInLoop
-
- enddo
-
- enddo kInLoop
-
-end subroutine InterpolateCubicSpline
-
-
-subroutine IntegrateCubicSpline(x,y,y2ndDer,n,x1,x2,y_integral)
-
-! Given the arrays x(1:n) and y(1:n), which tabulate a function,
-! and given the array y2ndDer(1:n), which is the output from
-! CubicSplineCoefficients above, this routine returns y_integral,
-! the integral of y from x1 to x2. The integration formula was
-! created by analytically integrating a cubic spline between each node.
-! This subroutine assumes that x is monotonically increasing, and
-! that x1 < x2.
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n ! number of nodes
- real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! location of nodes
- y, &! value at nodes
- y2ndDer ! dy^2/dx^2 at each node
- real(kind=RKIND), intent(in) :: &
- x1,x2 ! limits of integration
-
-! OUTPUT PARAMETERS:
-
- real(kind=RKIND), intent(out) :: &
- y_integral ! integral of y
-
-! local variables:
-
- integer :: i,j,k
- real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
-
- if (x1<x(1).or.x2>x(n).or.x1>x2) then
- print *, 'error on integration bounds'
- endif
-
- y_integral = 0.0
- eps1 = 1e-14*x2
-
- do j=1,n-1 ! loop through sections
- ! section x(j) ... x(j+1)
-
- if (x2<=x(j) +eps1) exit
- if (x1>=x(j+1)-eps1) cycle
-
- h = x(j+1) - x(j)
- h2 = h**2
-
- ! left side:
- if (x1<x(j)) then
- F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
- else
- A2 = (x(j+1)-x1 )**2/h2
- B2 = (x1 -x(j))**2/h2
- F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
- + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
- + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
- endif
-
- ! right side:
- if (x2>x(j+1)) then
- F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
- else
- A2 = (x(j+1)-x2 )**2/h2
- B2 = (x2 -x(j))**2/h2
- F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
- + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
- + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
- endif
-
- y_integral = y_integral + F2 - F1
-
- enddo ! j
-
- end subroutine IntegrateCubicSpline
-
-
- subroutine IntegrateColumnCubicSpline( &
- x,y,y2ndDer,n, &
- xOut,y_integral, nOut)
-
-! Given the arrays x(1:n) and y(1:n), which tabulate a function,
-! and given the array y2ndDer(1:n), which is the output from
-! CubicSplineCoefficients above, this routine returns
-! y_integral(1:nOut), the integral of y.
-! This is a cumulative integration, so that
-! y_integral(j) holds the integral of y from x(1) to xOut(j).
-! The integration formula was created by analytically integrating a
-! cubic spline between each node.
-! This subroutine assumes that both x and xOut are monotonically
-! increasing, and that all values of xOut are within the first and
-
-! INPUT PARAMETERS:
-
- integer, intent(in) :: &
- n, &! number of nodes
- nOut ! number of output locations to compute integral
- real(kind=RKIND), intent(in), dimension(n) :: &
- x, &! location of nodes
- y, &! value at nodes
- y2ndDer ! dy^2/dx^2 at each node
- real(kind=RKIND), dimension(nOut), intent(in) :: &
- xOut ! output locations to compute integral
-
-! OUTPUT PARAMETERS:
-
- real(kind=RKIND), dimension(nOut), intent(out) :: &
- y_integral ! integral from 0 to xOut
-
-! local variables:
-
- integer :: i,j,k
- real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
-
- y_integral = 0.0
- j = 1
- h = x(j+1) - x(j)
- h2 = h**2
- F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
- eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
-
- k_loop: do k = 1,nOut
-
- if (k>1) y_integral(k) = y_integral(k-1)
-
- do while(xOut(k) > x(j+1)-eps1)
- F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
-
- y_integral(k) = y_integral(k) + F2 - F1
- j = j+1
- h = x(j+1) - x(j)
- h2 = h**2
- F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
- if (abs(xOut(k) - x(j+1))<eps1) cycle k_loop
- enddo
-
- A2 = (x(j+1) - xOut(k))**2/h2
- B2 = (xOut(k) - x(j) )**2/h2
- F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
- + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
- + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
-
- y_integral(k) = y_integral(k) + F2 - F1
-
- if (k < nOut) then
- A2 = (x(j+1) -xOut(k))**2/h2
- B2 = (xOut(k) -x(j) )**2/h2
- F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
- + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
- + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
- endif
-
- enddo k_loop
-
- end subroutine IntegrateColumnCubicSpline
-
-
- subroutine InterpolateLinear( &
- x,y,n, &
- xOut,yOut,nOut)
-
-! Given the arrays x(1:n) and y(1:n), which tabulate a function,
-! this routine returns the linear interpolated values of yOut(1:nOut)
-! at xOut(1:nOut).
-! This subroutine assumes that both x and xOut are monotonically
-! increasing, and that all values of xOut are within the first and
-! last values of x.
-
-! !INPUT PARAMETERS:
-
- real (kind=RKIND), dimension(n), intent(in) :: &
- x, &! node location, input grid
- y ! interpolation variable, input grid
-
- real (kind=RKIND), dimension(nOut), intent(in) :: &
- xOut ! node location, output grid
-
- integer, intent(in) :: &
- N, &! number of nodes, input grid
- NOut ! number of nodes, output grid
-
-! !OUTPUT PARAMETERS:
-
- real (kind=RKIND), dimension(nOut), intent(out) :: &
- yOut ! interpolation variable, output grid
-
-!-----------------------------------------------------------------------
-!
-! local variables
-!
-!-----------------------------------------------------------------------
-
- integer :: &
- kIn, kOut ! counters
-
- kOut = 1
-
- kInLoop: do kIn = 1,n-1
-
- do while(xOut(kOut) < x(kIn+1))
-
- yOut(kOut) = y(kIn) &
- + (y(kIn+1)-y(kIn)) &
- /(x(kIn+1) -x(kIn) ) &
- *(xOut(kOut) -x(kIn) )
-
- kOut = kOut + 1
-
- if (kOut>nOut) exit kInLoop
-
- enddo
-
- enddo kInLoop
-
- end subroutine InterpolateLinear
-
-
- subroutine TestInterpolate
-
-! Test function to show how to operate the cubic spline subroutines
-
- integer, parameter :: &
- n = 10
- real (kind=RKIND), dimension(n) :: &
- y, x, y2ndDer
-
- integer, parameter :: &
- nOut = 100
- real (kind=RKIND), dimension(nOut) :: &
- yOut, xOut
-
- integer :: &
- k
-
-!-----------------------------------------------------------------------
-!
-! Create x, y, xOut
-!
-!-----------------------------------------------------------------------
-
- do k=1,n
- x(k) = k-4
- ! trig function:
- y(k) = sin(x(k)/2)
- enddo
-
- do k=1,nOut
- xOut(k) = x(1) + k/(nOut+1.0)*(x(n)-x(1))
- enddo
-
-!-----------------------------------------------------------------------
-!
-! Interpolate
-!
-!-----------------------------------------------------------------------
-
- ! First, compute second derivative values at each node, y2ndDer.
- call CubicSplineCoefficients(x,y,n,y2ndDer)
-
- ! Compute interpolated values yOut.
- call InterpolateCubicSpline( &
- x,y,y2ndDer,n, &
- xOut,yOut,nOut)
-
- ! The following output can be copied directly into Matlab
- print *, 'subplot(2,1,1)'
- print '(a,10f8.4,a)', 'x = [',x,'];'
- print '(a,10f8.4,a)', 'y = [',y,'];'
- print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
- print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
- print *, "plot(x,y,'-*r',xOut,yOut,'x')"
-
- ! Compute interpolated values yOut.
- call IntegrateColumnCubicSpline( &
- x,y,y2ndDer,n, &
- xOut,yOut,nOut)
-
- ! The following output can be copied directly into Matlab
- print *, 'subplot(2,1,2)'
- print '(a,10f8.4,a)', 'x = [',x,'];'
- print '(a,10f8.4,a)', 'y = 2*cos(-3/2) -2*cos(x/2);'
- print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
- print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
- print *, "plot(x,y,'-*r',xOut,yOut,'x')"
-
- end subroutine TestInterpolate
-
-end module spline_interpolation
-
Deleted: branches/source_renaming/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/source_renaming/src/operators/module_vector_reconstruction.F 2011-10-17 16:43:37 UTC (rev 1088)
+++ branches/source_renaming/src/operators/module_vector_reconstruction.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -1,196 +0,0 @@
-module vector_reconstruction
-
- use grid_types
- use configure
- use constants
- use RBF_interpolation
-
- implicit none
-
- public :: init_reconstruct, reconstruct
-
- contains
-
- subroutine init_reconstruct(grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: pre-compute coefficients used by the reconstruct() routine
- !
- ! Input: grid meta data
- !
- ! Output: grid % coeffs_reconstruct - coefficients used to reconstruct
- ! velocity vectors at cell centers
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(inout) :: grid
-
- ! temporary arrays needed in the (to be constructed) init procedure
- integer :: nCellsSolve
- integer, dimension(:,:), pointer :: edgesOnCell
- integer, dimension(:), pointer :: nEdgesOnCell
- integer :: i, iCell, iEdge, pointCount, maxEdgeCount
- real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
- real (kind=RKIND) :: r, cellCenter(3), alpha, tangentPlane(2,3)
- real (kind=RKIND), allocatable, dimension(:,:) :: edgeOnCellLocations, edgeOnCellNormals, &
- coeffs
-
- real(kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors
- real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
-
- real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-
- !========================================================
- ! arrays filled and saved during init procedure
- !========================================================
- coeffs_reconstruct => grid % coeffs_reconstruct % array
-
- !========================================================
- ! temporary variables needed for init procedure
- !========================================================
- xCell => grid % xCell % array
- yCell => grid % yCell % array
- zCell => grid % zCell % array
- xEdge => grid % xEdge % array
- yEdge => grid % yEdge % array
- zEdge => grid % zEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCellsSolve = grid % nCellsSolve
- edgeNormalVectors => grid % edgeNormalVectors % array
- cellTangentPlane => grid % cellTangentPlane % array
-
-
- ! init arrays
- coeffs_reconstruct = 0.0
-
- maxEdgeCount = maxval(nEdgesOnCell)
-
- allocate(edgeOnCellLocations(maxEdgeCount,3))
- allocate(edgeOnCellNormals(maxEdgeCount,3))
- allocate(coeffs(maxEdgeCount,3))
-
- ! loop over all cells to be solved on this block
- do iCell=1,nCellsSolve
- pointCount = nEdgesOnCell(iCell)
- cellCenter(1) = xCell(iCell)
- cellCenter(2) = yCell(iCell)
- cellCenter(3) = zCell(iCell)
-
- do i=1,pointCount
- iEdge = edgesOnCell(i,iCell)
- edgeOnCellLocations(i,1) = xEdge(iEdge)
- edgeOnCellLocations(i,2) = yEdge(iEdge)
- edgeOnCellLocations(i,3) = zEdge(iEdge)
- edgeOnCellNormals(i,:) = edgeNormalVectors(:, iEdge)
- end do
-
- alpha = 0.0
- do i=1,pointCount
- r = sqrt(sum((cellCenter - edgeOnCellLocations(i,:))**2))
- alpha = alpha + r
- enddo
- alpha = alpha/pointCount
-
- tangentPlane(1,:) = cellTangentPlane(:,1,iCell)
- tangentPlane(2,:) = cellTangentPlane(:,2,iCell)
-
- call rbfInterp_func_3DPlane_vec_const_dir_compCoeffs(pointCount, &
- edgeOnCellLocations(1:pointCount,:), &
- edgeOnCellNormals(1:pointCount,:), &
- cellCenter, alpha, tangentPlane, coeffs(1:pointCount,:))
-
- do i=1,pointCount
- coeffs_reconstruct(:,i,iCell) = coeffs(i,:)
- end do
-
- enddo ! iCell
-
- deallocate(edgeOnCellLocations)
- deallocate(edgeOnCellNormals)
- deallocate(coeffs)
-
- end subroutine init_reconstruct
-
- subroutine reconstruct(grid, u, uReconstructX, uReconstructY, uReconstructZ, uReconstructZonal, uReconstructMeridional)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: reconstruct vector field at cell centers based on radial basis functions
- !
- ! Input: grid meta data and vector component data residing at cell edges
- !
- ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND), dimension(:,:), intent(in) :: u
- real (kind=RKIND), dimension(:,:), intent(out) :: uReconstructX, uReconstructY, uReconstructZ
- real (kind=RKIND), dimension(:,:), intent(out) :: uReconstructZonal, uReconstructMeridional
-
- ! temporary arrays needed in the compute procedure
- integer :: nCellsSolve
- integer, dimension(:,:), pointer :: edgesOnCell
- integer, dimension(:), pointer :: nEdgesOnCell
- integer :: iCell,iEdge, i
- real(kind=RKIND), dimension(:), pointer :: latCell, lonCell
-
- real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
-
- logical :: on_a_sphere
-
- real (kind=RKIND) :: clat, slat, clon, slon
-
-
- ! stored arrays used during compute procedure
- coeffs_reconstruct => grid % coeffs_reconstruct % array
-
- ! temporary variables
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCellsSolve = grid % nCellsSolve
-
- latCell => grid % latCell % array
- lonCell => grid % lonCell % array
- on_a_sphere = grid % on_a_sphere
-
- ! init the intent(out)
- uReconstructX = 0.0
- uReconstructY = 0.0
- uReconstructZ = 0.0
-
- ! loop over cell centers
- do iCell=1,nCellsSolve
- ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
- ! in coeffs_reconstruct
- do i=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(i,iCell)
- uReconstructX(:,iCell) = uReconstructX(:,iCell) &
- + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
- uReconstructY(:,iCell) = uReconstructY(:,iCell) &
- + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
- uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
- + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
-
- enddo
- enddo ! iCell
-
- if(on_a_sphere) then
- do iCell=1,nCellsSolve
- clat = cos(latCell(iCell))
- slat = sin(latCell(iCell))
- clon = cos(lonCell(iCell))
- slon = sin(lonCell(iCell))
- uReconstructZonal(:,iCell) = -uReconstructX(:,iCell)*slon + uReconstructY(:,iCell)*clon
- uReconstructMeridional(:,iCell) = -(uReconstructX(:,iCell)*clon &
- + uReconstructY(:,iCell)*slon)*slat &
- + uReconstructZ(:,iCell)*clat
- end do
- else
- uReconstructZonal = uReconstructX
- uReconstructMeridional = uReconstructY
- end if
-
- end subroutine reconstruct
-
-end module vector_reconstruction
Copied: branches/source_renaming/src/operators/mpas_rbf_interpolation.F (from rev 1087, branches/source_renaming/src/operators/module_RBF_interpolation.F)
===================================================================
--- branches/source_renaming/src/operators/mpas_rbf_interpolation.F (rev 0)
+++ branches/source_renaming/src/operators/mpas_rbf_interpolation.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -0,0 +1,1824 @@
+module rbf_interpolation
+ use dmpar
+ use grid_types
+
+ implicit none
+ private
+ save
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Purpose: perform interpolation of scalar and vector functions in 2D
+! and 3D using Radial Basis Functions (RBFs).
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ ! Initialize the geometry that will be useful from interpolation
+ public :: mpas_rbf_interp_initialize
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Routines for perofrming interpolation in 2D (including Jacobian and Hessian)
+ ! at locations that vary using a function that is fixed. This is useful
+ ! 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 :: mpas_rbf_interp_loc_2D_sca_const_comp_coeffs, &
+ mpas_rbf_interp_loc_2D_sca_lin_comp_coeffs, &
+ mpas_rbf_interp_loc_2D_sca_const_eval_with_derivs, &
+ mpas_rbf_interp_loc_2D_sca_lin_eval_with_derivs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Routines for computing scalar interpolaiton coefficients in 3D (coplanar points
+ ! in 3D) with support for either constant or constant and linear basis
+ ! functions in addition to RBFs. In constrast to the two subroutines
+ ! 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.
+ ! 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(functionOrDerivAtSources*neumannCoefficients)
+ ! where functionOrDerivAtSource = functionAtSources where .not.(isInterface)
+ ! = functionNormalDerivAtSources where isInterface
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ public :: mpas_rbf_interp_func_3D_sca_const_dir_comp_coeffs, &
+ mpas_rbf_interp_func_3D_plane_sca_lin_dir_comp_coeffs, &
+ mpas_rbf_interp_func_3D_sca_lin_dir_comp_coeffs, &
+ mpas_rbf_interp_func_3D_sca_const_dir_neu_comp_coeffs, &
+ mpas_rbf_interp_func_3D_plane_sca_lin_dir_neu_comp_coeffs, &
+ mpas_rbf_interp_func_3D_sca_lin_dir_neu_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Routines for computing vector interpolaiton coefficients in 3D (coplanar points
+ ! in 3D) with support for only constant basis functions in addition to RBFs.
+ ! (Linear basis functions can cause problems: a linear vortex flow u = y xHat - x yHat
+ ! cannot be reconstructed from vectors that are only normal to the cell edges in 2D.
+ ! Therefore, we don't support them). As with the scalar 3D subroutines
+ ! 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.
+ ! The user supplies to these routines a set of sourcePoints and unitVectors
+ ! as well as a destinationPoint and, for the last 2 routines, flags for
+ ! which points are tangent to the interface and which of the supplied unitVectors
+ ! is the normal at the corresponding point.
+ !
+ ! The first two routines compute coefficients that can be used to interpolate
+ ! a vector function to the destination point given "function dot unitVector" values
+ ! at each source point. These routines are useful, for example, for reconstructing
+ ! the full vector velocity at cell centers from the normal components of the velocity
+ ! at cell faces (or cell edges in 2D), or for computing the full velocity at an
+ ! immersed boundary image point based on the normal velocity at several faces and
+ ! the full velocity at boundary points (e.g., a no-slip boundary condition).
+ !
+ ! The last two routines compute coefficients that can be used to interpolate
+ ! a vector function to the destination point given "function dot unitVector" values
+ ! at non-tangent source point and "dFunction/dn dot unitVector" values at
+ ! tangent source point. These routines are useful, for example, for computing the
+ ! full velocity at an immersed boundary image point based on the normal velocity at
+ ! several faces, the normal velocity at boundary points (e.g., u dot n = 0 for a
+ ! no-penetration boundary condition on a fixed boundary), and the normal derivative
+ ! of the tangential components of velocity at the boundary points (e.g., a free-slip
+ ! boundary condition).
+ ! Pseudocode for function reconstruction at the destinationPoint is as follows
+ ! dirichlet: functionAtDestination_i = sum_j(functionDotUnitVectorAtSources_j*coefficients_j,i)
+ ! for i = x,y,z
+ ! tangentNeumann: functionAtDestination_i &
+ ! = sum_(j where .not isTangent) (functionDotUnitVectorAtSources_j*coefficients_j,i) &
+ ! + sum_(j where isTangent) ((dFunction/dn dot UnitVector)_j*coefficients_j,i)
+ ! for i = x,y,z
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ public :: mpas_rbf_interp_func_3D_vec_const_dir_comp_coeffs, &
+ mpas_rbf_interp_func_3D_plane_vec_const_dir_comp_coeffs!, &
+ !mpas_rbf_interp_func_3d_vec_const_tan_neu_comp_coeffs, &
+ !mpas_rbf_interp_func_3d_plane_vec_const_tan_neu_comp_coeffs
+
+ contains
+
+ subroutine mpas_rbf_interp_initialize(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: compute geometric fields that will be potentially useful for calling
+ ! the interpolation routines
+ !
+ ! Input: the grid
+ !
+ ! Output:
+ ! edgeNormalVectors - the unit vector at the center of each edge tangent to the sphere
+ ! cellTangentPlane - 2 orthogonal unit vectors in the tangent plane of each cell
+ ! The first unit vector is chosen to point toward the center of the first
+ ! edge on the cell.
+ ! localVerticalUnitVectors - the unit normal vector of the tangent plane at the center
+ ! of each cell
+ !
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+
+ integer :: nCells, nEdges
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
+ integer :: iEdge, iCell, cell1, cell2
+ real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
+ real(kind=RKIND), dimension(:,:), pointer :: localVerticalUnitVectors, edgeNormalVectors
+ real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
+ real(kind=RKIND), dimension(3) :: xHatPlane, yHatPlane, rHat
+ real(kind=RKIND) :: normalDotRHat
+ logical :: on_a_sphere
+
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ xEdge => grid % xEdge % array
+ yEdge => grid % yEdge % array
+ zEdge => grid % zEdge % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ nCells = grid % nCells
+ nEdges = grid % nEdges
+ on_a_sphere = grid % on_a_sphere
+
+ localVerticalUnitVectors => grid % localVerticalUnitVectors % array
+ edgeNormalVectors => grid % edgeNormalVectors % array
+ cellTangentPlane => grid % cellTangentPlane % array
+
+ ! init arrays
+ edgeNormalVectors = 0
+ localVerticalUnitVectors = 0
+
+ ! loop over all cells to be solved on this block
+ do iCell=1,nCells
+ if(on_a_sphere) then
+ localVerticalUnitVectors(1,iCell) = xCell(iCell)
+ localVerticalUnitVectors(2,iCell) = yCell(iCell)
+ localVerticalUnitVectors(3,iCell) = zCell(iCell)
+ call mpas_unit_vec_in_r3(localVerticalUnitVectors(:,iCell))
+ else ! on a plane
+ localVerticalUnitVectors(:,iCell) = (/ 0., 0., 1. /)
+ end if
+ end do
+
+ do iEdge = 1,nEdges
+ iCell = cellsOnEdge(1,iEdge) ! the normal vector points from the first cell toward the edge
+ if(iCell == nCells+1) then ! this is a boundary edge
+ ! the first cell bordering this edge is not real, use the second cell
+ ! The normal should always point outward at boundaries, away from the valid cell center
+ iCell = cellsOnEdge(2,iEdge)
+ end if
+ ! the normal points from the cell location to the edge location
+ edgeNormalVectors(1,iEdge) = xEdge(iEdge) - xCell(iCell)
+ edgeNormalVectors(2,iEdge) = yEdge(iEdge) - yCell(iCell)
+ edgeNormalVectors(3,iEdge) = zEdge(iEdge) - zCell(iCell)
+ call mpas_unit_vec_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
+ rHat = localVerticalUnitVectors(:,iCell)
+ normalDotRHat = sum(edgeNormalVectors(:,iEdge)*rHat)
+ xHatPlane = edgeNormalVectors(:,iEdge) - normalDotRHat*rHat
+ call mpas_unit_vec_in_r3(xHatPlane)
+
+ call mpas_cross_product_in_r3(rHat, xHatPlane, yHatPlane)
+ call mpas_unit_vec_in_r3(yHatPlane) ! just to be sure...
+ cellTangentPlane(:,1,iCell) = xHatPlane
+ cellTangentPlane(:,2,iCell) = yHatPlane
+ end do
+
+ end subroutine mpas_rbf_interp_initialize
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 2D that can be used to
+ ! reconstruct a given scalar function at varying locations. This is useful
+ ! 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.
+ ! The reconstruction is performed with basis functions that are RBFs and constant
+ ! Input:
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! coeffCount - the size of coefficients, must be at least pointCount + 1
+ ! points - the location of the "source" points in the 2D space where the values of
+ ! the function are known
+ ! fieldValues - the values of the function of interest at the points
+ ! 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 needed to perform interpolation of the funciton
+ ! at destination points yet to be specified
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mpas_rbf_interp_loc_2d_sca_const_comp_coeffs(pointCount, coeffCount, &
+ points, fieldValues, alpha, coefficients)
+
+ integer, intent(in) :: pointCount, coeffCount
+ real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
+ real(kind=RKIND), dimension(pointCount), intent(in) :: fieldValues
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(coeffCount), intent(out) :: coefficients
+
+ integer :: i, j, matrixSize
+ real(kind=RKIND), dimension(pointCount+1,pointCount+1) :: matrix
+ real(kind=RKIND), dimension(pointCount+1) :: rhs
+ integer, dimension(pointCount+1) :: pivotIndices
+ real(kind=RKIND) :: rSquared
+
+ matrixSize = pointCount+1
+ coefficients = 0.0
+ matrix = 0.0
+ rhs = 0.0
+
+ 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
+ do j=1,pointCount
+ matrix(pointCount+1,j) = 1.0
+ matrix(j,pointCount+1) = 1.0
+ end do
+
+ call mpas_legs(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
+ coefficients(1:matrixSize), pivotIndices(1:matrixSize))
+
+ end subroutine mpas_rbf_interp_loc_2d_sca_const_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Compute interpolation coefficients in 2D that can be used to
+ ! reconstruct a given scalar function at varying locations. This is useful
+ ! 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.
+ ! The reconstruction is performed with basis functions that are RBFs plus constant
+ ! and linear
+ ! Input:
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! coeffCount - the size of coefficients, must be at least pointCount + 3
+ ! points - the location of the "source" points in the 2D space where the values of
+ ! the function are known
+ ! fieldValues - the values of the function of interest at the points
+ ! 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 needed to perform interpolation of the funciton
+ ! at destination points yet to be specified
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mpas_rbf_interp_loc_2d_sca_lin_comp_coeffs(pointCount, coeffCount, &
+ points, fieldValues, alpha, coefficients)
+
+ integer, intent(in) :: pointCount, coeffCount
+ real(kind=RKIND), dimension(pointCount,2), intent(in) :: points
+ real(kind=RKIND), dimension(pointCount), intent(in) :: fieldValues
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(coeffCount), intent(out) :: coefficients
+
+ integer :: i, j, matrixSize
+ real(kind=RKIND), dimension(pointCount+3,pointCount+3) :: matrix
+ real(kind=RKIND), dimension(pointCount+3) :: rhs
+ integer, dimension(pointCount+3) :: pivotIndices
+ real(kind=RKIND) :: rSquared
+
+ coefficients = 0.0
+ matrix = 0.0
+ rhs = 0.0
+
+ 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
+ matrixSize = pointCount+3
+ do j=1,pointCount
+ matrix(pointCount+1,j) = 1.0
+ 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
+ call mpas_legs(matrix(1:matrixSize,1:matrixSize), matrixSize, rhs(1:matrixSize), &
+ coefficients(1:matrixSize), pivotIndices(1:matrixSize))
+
+ end subroutine mpas_rbf_interp_loc_2d_sca_lin_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Evalute a scalar function in 2D using coefficients computed in
+ ! 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
+ ! RBF reconstruction of a function (e.g., a height field) that minimizes
+ ! the distance to a point in 3D space. The reconstruction is performed
+ ! with basis functions that are RBFs and constant
+ ! Input:
+ ! fieldCount - the number fields to be evaluated. This is useful for reconstructing,
+ ! for example, the x-, y- and z-components of a vector field at the same
+ ! point in 2D
+ ! coeffCount - the size of coefficients, must be at least pointCount + 1
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! coefficients - the coefficients needed to perform interpolation of the funciton
+ ! at the evaluationPoint
+ ! evaluationPoint - the point in 2D where the function is to be reconstructed
+ ! points - the location of the "source" points in the 2D space where the values of
+ ! the function are known
+ ! alpha - a constant that give the characteristic length scale of the RBFs,
+ ! should be on the order of the distance between points
+ ! Output:
+ ! derivs - the value of the function, the 2 components of its Jacobian and
+ ! the 3 unique components of its Hessian at the evaluationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mpas_rbf_interp_loc_2d_sca_const_eval_with_derivs(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
+ real(kind=RKIND), intent(in) :: alpha
+
+ real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
+
+ integer :: pointIndex
+ real(kind=RKIND) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+ derivs = 0.0
+ do pointIndex = 1, pointCount
+ x = (evaluationPoint(1) - points(pointIndex,1))
+ y = (evaluationPoint(2) - points(pointIndex,2))
+ rSquared = x**2 + y**2
+ call mpas_evaluate_rbf_and_derivs(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 mpas_evaluate_rbf_and_derivs(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
+ derivs(1,:) = derivs(1,:) + coefficients(pointCount+1,:)
+ end subroutine mpas_rbf_interp_loc_2d_sca_const_eval_with_derivs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: Evalute a scalar function in 2D using coefficients computed in
+ ! 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
+ ! RBF reconstruction of a function (e.g., a height field) that minimizes
+ ! the distance to a point in 3D space. The reconstruction is performed
+ ! with basis functions that are RBFs, constant and linear
+ ! Input:
+ ! fieldCount - the number fields to be evaluated. This is useful for reconstructing,
+ ! for example, the x-, y- and z-components of a vector field at the same
+ ! point in 2D
+ ! coeffCount - the size of coefficients, must be at least pointCount + 1
+ ! pointCount - the number of "source" points and functionValues supplied
+ ! coefficients - the coefficients needed to perform interpolation of the funciton
+ ! at the evaluationPoint
+ ! evaluationPoint - the point in 2D where the function is to be reconstructed
+ ! points - the location of the "source" points in the 2D space where the values of
+ ! the function are known
+ ! alpha - a constant that give the characteristic length scale of the RBFs,
+ ! should be on the order of the distance between points
+ ! Output:
+ ! derivs - the value of the function, the 2 components of its Jacobian and
+ ! the 3 unique components of its Hessian at the evaluationPoint
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine mpas_rbf_interp_loc_2d_sca_lin_eval_with_derivs(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
+ real(kind=RKIND), intent(in) :: alpha
+
+ real(kind=RKIND), dimension(6,fieldCount), intent(out) :: derivs
+
+ integer :: pointIndex
+ real(kind=RKIND) :: x, y, rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+ derivs = 0.0
+ do pointIndex = 1, pointCount
+ x = (evaluationPoint(1) - points(pointIndex,1))
+ y = (evaluationPoint(2) - points(pointIndex,2))
+ rSquared = x**2 + y**2
+ call mpas_evaluate_rbf_and_derivs(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 mpas_evaluate_rbf_and_derivs(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
+ 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,:)
+
+ end subroutine mpas_rbf_interp_loc_2d_sca_lin_eval_with_derivs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_sca_const_dir_comp_coeffs( &
+ 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 mpas_set_up_scalar_rbf_dirichlet_matrix_and_rhs(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 mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_sca_const_dir_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_plane_sca_lin_dir_comp_coeffs( &
+ 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 mpas_set_up_scalar_rbf_dirichlet_matrix_and_rhs(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 mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_plane_sca_lin_dir_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_sca_lin_dir_comp_coeffs(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 mpas_set_up_scalar_rbf_dirichlet_matrix_and_rhs(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 mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ coefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_sca_lin_dir_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_sca_const_dir_neu_comp_coeffs( &
+ pointCount, sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+1 !! 1 extra space for constant
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(neumannMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(rhsCopy(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ neumannMatrix = 0.0
+ rhs = 0.0
+ rhsCopy = 0.0
+ coeffs = 0.0
+
+ call mpas_set_up_scalar_rbf_matrix_and_rhs(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
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0
+ 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
+
+ ! solve each linear system
+ rhsCopy = rhs
+ call mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call mpas_legs(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(neumannMatrix)
+ deallocate(rhs)
+ deallocate(rhsCopy)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_sca_const_dir_neu_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_plane_sca_lin_dir_neu_comp_coeffs( &
+ pointCount, sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, planeBasisVectors, dirichletCoefficients, neumannCoefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+ 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) :: &
+ dirichletCoefficients, neumannCoefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+3 !! 3 extra space for constant and 2 planar dimensions
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(neumannMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(rhsCopy(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ neumannMatrix = 0.0
+ rhs = 0.0
+ rhsCopy = 0.0
+ coeffs = 0.0
+
+ call mpas_set_up_scalar_rbf_matrix_and_rhs(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
+ dirichletMatrix(i,pointCount+2) = sum(sourcePoints(i,1:3)*planeBasisVectors(1,:))
+ dirichletMatrix(i,pointCount+3) = sum(sourcePoints(i,1:3)*planeBasisVectors(2,:))
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0
+ neumannMatrix(i,pointCount+2) = sum(interfaceNormals(i,1:3)*planeBasisVectors(1,:))
+ neumannMatrix(i,pointCount+3) = sum(interfaceNormals(i,1:3)*planeBasisVectors(2,:))
+ else
+ neumannMatrix(i,pointCount+1:pointCount+3) &
+ = dirichletMatrix(i,pointCount+1:pointCount+3)
+ end if
+ dirichletMatrix(pointCount+1:pointCount+3,i) &
+ = dirichletMatrix(i,pointCount+1:pointCount+3)
+ neumannMatrix(pointCount+1:pointCount+3,i) &
+ = neumannMatrix(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
+ rhsCopy = rhs
+ call mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call mpas_legs(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(neumannMatrix)
+ deallocate(rhs)
+ deallocate(rhsCopy)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_plane_sca_lin_dir_neu_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_sca_lin_dir_neu_comp_coeffs(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletCoefficients, neumannCoefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount), intent(out) :: &
+ dirichletCoefficients, neumannCoefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND), dimension(:,:), pointer :: dirichletMatrix, neumannMatrix
+ real(kind=RKIND), dimension(:), pointer :: rhs, rhsCopy, coeffs
+ integer, dimension(:), pointer :: pivotIndices
+
+ matrixSize = pointCount+4 !! 4 extra space for constant and linear in 3D
+
+ allocate(dirichletMatrix(matrixSize,matrixSize))
+ allocate(neumannMatrix(matrixSize,matrixSize))
+ allocate(rhs(matrixSize))
+ allocate(rhsCopy(matrixSize))
+ allocate(coeffs(matrixSize))
+ allocate(pivotIndices(matrixSize))
+
+ dirichletMatrix = 0.0
+ neumannMatrix = 0.0
+ rhs = 0.0
+ rhsCopy = 0.0
+ coeffs = 0.0
+
+ call mpas_set_up_scalar_rbf_matrix_and_rhs(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
+ dirichletMatrix(i,pointCount+2:pointCount+4) = sourcePoints(i,1:3)
+ if(isInterface(i)) then
+ neumannMatrix(i,pointCount+1) = 0.0
+ 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
+ rhs(pointCount+2:pointCount+4) = destinationPoint(1:3)
+
+ ! solve each linear system
+ rhsCopy = rhs
+ call mpas_legs(dirichletMatrix, matrixSize, rhs, coeffs, pivotIndices)
+ dirichletCoefficients = coeffs(1:pointCount)
+
+ call mpas_legs(neumannMatrix, matrixSize, rhsCopy, coeffs, pivotIndices)
+ neumannCoefficients = coeffs(1:pointCount)
+
+ deallocate(dirichletMatrix)
+ deallocate(neumannMatrix)
+ deallocate(rhs)
+ deallocate(rhsCopy)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_sca_lin_dir_neu_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_vec_const_dir_comp_coeffs(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+
+ 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 mpas_set_up_vector_dirichlet_rbf_matrix_and_rhs(pointCount, 3, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+
+ do i = 1, pointCount
+ matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+ matrix(pointCount+1:pointCount+3,i) &
+ = matrix(i,pointCount+1:pointCount+3)
+ end do
+ do i = 1, 3
+ rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+ end do
+
+ ! solve each linear system
+ do i = 1, 3
+ matrixCopy = matrix
+ call mpas_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 mpas_rbf_interp_func_3d_vec_const_dir_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_plane_vec_const_dir_comp_coeffs(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, planeBasisVectors, coefficients)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+ real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+
+ integer :: i, j
+ integer :: matrixSize
+
+ real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
+
+ 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 mpas_set_up_vector_dirichlet_rbf_matrix_and_rhs(pointCount, 2, &
+ planarSourcePoints, planarUnitVectors, planarDestinationPoint, &
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+
+ do i = 1, pointCount
+ matrix(i,pointCount+1:pointCount+2) = planarUnitVectors(i,:)
+ 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 mpas_legs(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+ call mpas_legs(matrixCopy, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+
+
+ do i = 1,3
+ coefficients(:,i) = planeBasisVectors(1,i)*coeffs(1:pointCount,1) &
+ + planeBasisVectors(2,i)*coeffs(1:pointCount,2)
+ end do
+
+ deallocate(matrix)
+ deallocate(matrixCopy)
+ deallocate(rhs)
+ deallocate(coeffs)
+ deallocate(pivotIndices)
+
+ end subroutine mpas_rbf_interp_func_3d_plane_vec_const_dir_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_vec_const_tan_neu_comp_coeffs(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 mpas_set_up_vector_free_slip_rbf_matrix_and_rhs(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 mpas_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 mpas_rbf_interp_func_3d_vec_const_tan_neu_comp_coeffs
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 mpas_rbf_interp_func_3d_plane_vec_const_tan_neu_comp_coeffs(&
+ 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 mpas_set_up_vector_free_slip_rbf_matrix_and_rhs(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 mpas_legs(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+ call mpas_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 mpas_rbf_interp_func_3d_plane_vec_const_tan_neu_comp_coeffs
+
+
+!!!!!!!!!!!!!!!!!!!!!
+! private subroutines
+!!!!!!!!!!!!!!!!!!!!!
+
+ function evaluateRBF(rSquared) result(rbfValue)
+ real(kind=RKIND), intent(in) :: rSquared
+ real(kind=RKIND) :: rbfValue
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+
+ end function evaluateRBF
+
+ subroutine mpas_evaluate_rbf_and_deriv(rSquared, rbfValue, rbfDerivOverR)
+ real(kind=RKIND), intent(in) :: rSquared
+ real(kind=RKIND), intent(out) :: rbfValue, rbfDerivOverR
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+ rbfDerivOverR = -rbfValue**3
+
+ end subroutine mpas_evaluate_rbf_and_deriv
+
+ subroutine mpas_evaluate_rbf_and_derivs(rSquared, rbfValue, rbfDerivOverR, rbfSecondDeriv)
+ real(kind=RKIND), intent(in) :: rSquared
+ real(kind=RKIND), intent(out) :: rbfValue, rbfDerivOverR, rbfSecondDeriv
+
+ ! inverse multiquadratic
+ rbfValue = 1/sqrt(1 + rSquared)
+ rbfDerivOverR = -rbfValue**3
+ rbfSecondDeriv = (2*rSquared-1)*rbfValue**5
+
+ end subroutine mpas_evaluate_rbf_and_derivs
+
+ subroutine mpas_set_up_scalar_rbf_dirichlet_matrix_and_rhs(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 mpas_set_up_scalar_rbf_dirichlet_matrix_and_rhs
+
+ subroutine mpas_set_up_scalar_rbf_matrix_and_rhs(pointCount, &
+ sourcePoints, isInterface, interfaceNormals, destinationPoint, &
+ alpha, dirichletMatrix, neumannMatrix, rhs)
+
+ integer, intent(in) :: pointCount
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isInterface
+ real(kind=RKIND), dimension(pointCount,3), intent(in) :: interfaceNormals
+ real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: &
+ dirichletMatrix, neumannMatrix
+ real(kind=RKIND), dimension(pointCount), intent(out) :: rhs
+
+ integer :: i, j
+
+ real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalDotX
+
+ do j = 1, pointCount
+ if(isInterface(j)) then
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalDotX = sum(interfaceNormals(j,:) &
+ * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call mpas_evaluate_rbf_and_deriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ dirichletMatrix(i,j) = rbfValue
+ neumannMatrix(i,j) = rbfDerivOverR*normalDotX
+ end do
+ else
+ 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 if
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j) = evaluateRBF(rSquared)
+ end do
+
+ end subroutine mpas_set_up_scalar_rbf_matrix_and_rhs
+
+ subroutine mpas_set_up_vector_dirichlet_rbf_matrix_and_rhs(pointCount, dimensions, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
+
+ integer, intent(in) :: pointCount, dimensions
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
+
+ integer :: i, j
+
+ real(kind=RKIND) :: rSquared, rbfValue, unitVectorDotProduct
+
+ do j = 1, pointCount
+ do i = j, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfValue*unitVectorDotProduct
+ matrix(j,i) = matrix(i,j)
+ end do
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
+
+ end subroutine mpas_set_up_vector_dirichlet_rbf_matrix_and_rhs
+
+ subroutine mpas_set_up_vector_free_slip_rbf_matrix_and_rhs(pointCount, dimensions, &
+ sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+ alpha, matrix, rhs)
+
+ integer, intent(in) :: pointCount, dimensions
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: sourcePoints
+ logical, dimension(pointCount), intent(in) :: isTangentToInterface
+ integer, dimension(pointCount), intent(in) :: normalVectorIndex
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(in) :: unitVectors
+ real(kind=RKIND), dimension(dimensions), intent(in) :: destinationPoint
+ real(kind=RKIND), intent(in) :: alpha
+ real(kind=RKIND), dimension(pointCount,pointCount), intent(out) :: matrix
+ real(kind=RKIND), dimension(pointCount,dimensions), intent(out) :: rhs
+
+ integer :: i, j
+
+ real(kind=RKIND) :: rSquared, rbfValue, rbfDerivOverR, normalVector(dimensions), &
+ normalDotX, unitVectorDotProduct
+
+ do j = 1, pointCount
+ if(isTangentToInterface(j)) then
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ normalVector = unitVectors(normalVectorIndex(j),:)
+ normalDotX = sum(normalVector * (sourcePoints(j,:)-sourcePoints(i,:)))
+ call mpas_evaluate_rbf_and_deriv(rSquared, rbfValue, rbfDerivOverR)
+ rbfDerivOverR = rbfDerivOverR/alpha**2
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfDerivOverR*normalDotX*unitVectorDotProduct
+ end do
+ else
+ do i = 1, pointCount
+ rSquared = sum((sourcePoints(i,:)-sourcePoints(j,:))**2)/alpha**2
+ rbfValue = evaluateRBF(rSquared)
+ unitVectorDotProduct = sum(unitVectors(i,:)*unitVectors(j,:))
+ matrix(i,j) = rbfValue*unitVectorDotProduct
+ end do
+ end if
+ end do
+
+ do j = 1, pointCount
+ rSquared = sum((destinationPoint-sourcePoints(j,:))**2)/alpha**2
+ rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
+ end do
+
+ end subroutine mpas_set_up_vector_free_slip_rbf_matrix_and_rhs
+
+ subroutine mpas_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 mpas_unit_vec_in_r3
+
+ subroutine mpas_cross_product_in_r3(p_1,p_2,p_out)
+ real (kind=RKIND), intent(in) :: p_1 (3), p_2 (3)
+ 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 mpas_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 mpas_legs (A,N,B,X,INDX)
+!
+! WRITE (6, "(F16.8)") (X(I), I=1,N)
+!END PROGRAM EX43
+
+
+subroutine mpas_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, INTENT (IN) :: N
+ integer :: I,J
+ integer, INTENT (OUT), DIMENSION (N) :: INDX
+ real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+ real(kind=RKIND), INTENT (INOUT), DIMENSION (N) :: B
+ real(kind=RKIND), 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 mpas_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, INTENT (IN) :: N
+ integer :: I,J,K
+ integer, INTENT (OUT), DIMENSION (N) :: INDX
+ real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N):: A
+ real(kind=RKIND), INTENT (OUT), DIMENSION (N,N):: X
+ real(kind=RKIND), 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, INTENT (IN) :: N
+ integer :: I,J,K,ITMP
+ integer, INTENT (OUT), DIMENSION (N) :: INDX
+ real(kind=RKIND) :: C1,PI,PI1,PJ
+ real(kind=RKIND), INTENT (INOUT), DIMENSION (N,N) :: A
+ real(kind=RKIND), 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 rbf_interpolation
+
Copied: branches/source_renaming/src/operators/mpas_spline_interpolation.F (from rev 1087, branches/source_renaming/src/operators/module_spline_interpolation.F)
===================================================================
--- branches/source_renaming/src/operators/mpas_spline_interpolation.F (rev 0)
+++ branches/source_renaming/src/operators/mpas_spline_interpolation.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -0,0 +1,430 @@
+module spline_interpolation
+
+ implicit none
+
+ private
+
+ public :: mpas_cubic_spline_coefficients, &
+ mpas_interpolate_cubic_spline, &
+ mpas_integrate_cubic_spline, &
+ mpas_integrate_column_cubic_spline, &
+ mpas_interpolate_linear, &
+ mpas_test_interpolate
+
+! Short Descriptions:
+
+! mpas_cubic_spline_coefficients: Compute second derivatives at nodes.
+! This must be run before any of the other cubic spine functions.
+
+! mpas_interpolate_cubic_spline: Compute cubic spline interpolation.
+
+! mpas_integrate_cubic_spline: Compute a single integral from spline data.
+
+! mpas_integrate_column_cubic_spline: Compute multiple integrals from spline data.
+
+! mpas_interpolate_linear: Compute linear interpolation.
+
+! mpas_test_interpolate: Test spline interpolation subroutines.
+
+ contains
+
+ subroutine mpas_cubic_spline_coefficients(x,y,n,y2ndDer)
+
+! Given arrays x(1:n) and y(1:n) containing a function,
+! i.e., y(i) = f(x(i)), with x monotonically increasing
+! this routine returns an array y2ndDer(1:n) that contains
+! the second derivatives of the interpolating function at x(1:n).
+! This routine uses boundary conditions for a natural spline,
+! with zero second derivative on that boundary.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y ! value at nodes
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out), dimension(n) :: &
+ y2ndDer ! dy^2/dx^2 at each node
+
+! local variables:
+
+ integer :: i
+ real(kind=RKIND) :: &
+ temp,xRatio,a(n)
+
+ y2ndDer(1)=0.0
+ y2ndDer(n)=0.0
+ a(1)=0.0
+
+ do i=2,n-1
+ xRatio=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ temp=1.0/(2.0+xRatio*y2ndDer(i-1))
+ y2ndDer(i)=temp*(xRatio-1.0)
+ a(i) = temp*(6.0*((y(i+1)-y(i))/(x(i+1)-x(i)) &
+ -(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1)) &
+ -xRatio*a(i-1))
+ enddo
+
+ do i=n-1,1,-1
+ y2ndDer(i)=y2ndDer(i)*y2ndDer(i+1)+a(i)
+ enddo
+
+ end subroutine mpas_cubic_spline_coefficients
+
+
+ subroutine mpas_interpolate_cubic_spline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns the
+! cubic-spline interpolated values of yOut(1:nOut) at xOut(1:nOut).
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+! last values of x.
+
+! INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y, &! interpolation variable, input grid
+ y2ndDer ! 2nd derivative of y at nodes
+
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
+
+ integer, intent(in) :: &
+ n, &! number of nodes, input grid
+ nOut ! number of nodes, output grid
+
+! OUTPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
+
+! local variables:
+
+ integer :: &
+ kIn, kOut ! counters
+
+ real (kind=RKIND) :: &
+ a, b, h
+
+ kOut = 1
+
+ kInLoop: do kIn = 1,n-1
+
+ h = x(kIn+1)-x(kIn)
+
+ do while(xOut(kOut) < x(kIn+1))
+
+ a = (x(kIn+1)-xOut(kOut))/h
+ b = (xOut(kOut)-x (kIn) )/h
+ yOut(kOut) = a*y(kIn) + b*y(kIn+1) &
+ + ((a**3-a)*y2ndDer(kIn) + (b**3-b)*y2ndDer(kIn+1)) &
+ *(h**2)/6.0
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+end subroutine mpas_interpolate_cubic_spline
+
+
+subroutine mpas_integrate_cubic_spline(x,y,y2ndDer,n,x1,x2,y_integral)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns y_integral,
+! the integral of y from x1 to x2. The integration formula was
+! created by analytically integrating a cubic spline between each node.
+! This subroutine assumes that x is monotonically increasing, and
+! that x1 < x2.
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n ! number of nodes
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y, &! value at nodes
+ y2ndDer ! dy^2/dx^2 at each node
+ real(kind=RKIND), intent(in) :: &
+ x1,x2 ! limits of integration
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), intent(out) :: &
+ y_integral ! integral of y
+
+! local variables:
+
+ integer :: i,j,k
+ real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
+
+ if (x1<x(1).or.x2>x(n).or.x1>x2) then
+ print *, 'error on integration bounds'
+ endif
+
+ y_integral = 0.0
+ eps1 = 1e-14*x2
+
+ do j=1,n-1 ! loop through sections
+ ! section x(j) ... x(j+1)
+
+ if (x2<=x(j) +eps1) exit
+ if (x1>=x(j+1)-eps1) cycle
+
+ h = x(j+1) - x(j)
+ h2 = h**2
+
+ ! left side:
+ if (x1<x(j)) then
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ else
+ A2 = (x(j+1)-x1 )**2/h2
+ B2 = (x1 -x(j))**2/h2
+ F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ ! right side:
+ if (x2>x(j+1)) then
+ F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
+ else
+ A2 = (x(j+1)-x2 )**2/h2
+ B2 = (x2 -x(j))**2/h2
+ F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ y_integral = y_integral + F2 - F1
+
+ enddo ! j
+
+ end subroutine mpas_integrate_cubic_spline
+
+
+ subroutine mpas_integrate_column_cubic_spline( &
+ x,y,y2ndDer,n, &
+ xOut,y_integral, nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! and given the array y2ndDer(1:n), which is the output from
+! CubicSplineCoefficients above, this routine returns
+! y_integral(1:nOut), the integral of y.
+! This is a cumulative integration, so that
+! y_integral(j) holds the integral of y from x(1) to xOut(j).
+! The integration formula was created by analytically integrating a
+! cubic spline between each node.
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+
+! INPUT PARAMETERS:
+
+ integer, intent(in) :: &
+ n, &! number of nodes
+ nOut ! number of output locations to compute integral
+ real(kind=RKIND), intent(in), dimension(n) :: &
+ x, &! location of nodes
+ y, &! value at nodes
+ y2ndDer ! dy^2/dx^2 at each node
+ real(kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! output locations to compute integral
+
+! OUTPUT PARAMETERS:
+
+ real(kind=RKIND), dimension(nOut), intent(out) :: &
+ y_integral ! integral from 0 to xOut
+
+! local variables:
+
+ integer :: i,j,k
+ real(kind=RKIND) :: h,h2, A2,B2, F1,F2, eps1
+
+ y_integral = 0.0
+ j = 1
+ h = x(j+1) - x(j)
+ h2 = h**2
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ eps1 = 0.0 ! note: could use 1e-12*xOut(nOut)
+
+ k_loop: do k = 1,nOut
+
+ if (k>1) y_integral(k) = y_integral(k-1)
+
+ do while(xOut(k) > x(j+1)-eps1)
+ F2 = y(j+1)*h*0.5 - y2ndDer(j+1)*h**3/24.0
+
+ y_integral(k) = y_integral(k) + F2 - F1
+ j = j+1
+ h = x(j+1) - x(j)
+ h2 = h**2
+ F1 = -y(j)*h*0.5 + y2ndDer(j)*h**3/24.0
+ if (abs(xOut(k) - x(j+1))<eps1) cycle k_loop
+ enddo
+
+ A2 = (x(j+1) - xOut(k))**2/h2
+ B2 = (xOut(k) - x(j) )**2/h2
+ F2 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+
+ y_integral(k) = y_integral(k) + F2 - F1
+
+ if (k < nOut) then
+ A2 = (x(j+1) -xOut(k))**2/h2
+ B2 = (xOut(k) -x(j) )**2/h2
+ F1 = 0.5*h*( -y(j)*A2 + y(j+1)*B2 &
+ + y2ndDer(j) *h2*(-0.5*A2**2 + A2)/6.0 &
+ + y2ndDer(j+1)*h2*( 0.5*B2**2 - B2)/6.0 )
+ endif
+
+ enddo k_loop
+
+ end subroutine mpas_integrate_column_cubic_spline
+
+
+ subroutine mpas_interpolate_linear( &
+ x,y,n, &
+ xOut,yOut,nOut)
+
+! Given the arrays x(1:n) and y(1:n), which tabulate a function,
+! this routine returns the linear interpolated values of yOut(1:nOut)
+! at xOut(1:nOut).
+! This subroutine assumes that both x and xOut are monotonically
+! increasing, and that all values of xOut are within the first and
+! last values of x.
+
+! !INPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(n), intent(in) :: &
+ x, &! node location, input grid
+ y ! interpolation variable, input grid
+
+ real (kind=RKIND), dimension(nOut), intent(in) :: &
+ xOut ! node location, output grid
+
+ integer, intent(in) :: &
+ N, &! number of nodes, input grid
+ NOut ! number of nodes, output grid
+
+! !OUTPUT PARAMETERS:
+
+ real (kind=RKIND), dimension(nOut), intent(out) :: &
+ yOut ! interpolation variable, output grid
+
+!-----------------------------------------------------------------------
+!
+! local variables
+!
+!-----------------------------------------------------------------------
+
+ integer :: &
+ kIn, kOut ! counters
+
+ kOut = 1
+
+ kInLoop: do kIn = 1,n-1
+
+ do while(xOut(kOut) < x(kIn+1))
+
+ yOut(kOut) = y(kIn) &
+ + (y(kIn+1)-y(kIn)) &
+ /(x(kIn+1) -x(kIn) ) &
+ *(xOut(kOut) -x(kIn) )
+
+ kOut = kOut + 1
+
+ if (kOut>nOut) exit kInLoop
+
+ enddo
+
+ enddo kInLoop
+
+ end subroutine mpas_interpolate_linear
+
+
+ subroutine mpas_test_interpolate
+
+! Test function to show how to operate the cubic spline subroutines
+
+ integer, parameter :: &
+ n = 10
+ real (kind=RKIND), dimension(n) :: &
+ y, x, y2ndDer
+
+ integer, parameter :: &
+ nOut = 100
+ real (kind=RKIND), dimension(nOut) :: &
+ yOut, xOut
+
+ integer :: &
+ k
+
+!-----------------------------------------------------------------------
+!
+! Create x, y, xOut
+!
+!-----------------------------------------------------------------------
+
+ do k=1,n
+ x(k) = k-4
+ ! trig function:
+ y(k) = sin(x(k)/2)
+ enddo
+
+ do k=1,nOut
+ xOut(k) = x(1) + k/(nOut+1.0)*(x(n)-x(1))
+ enddo
+
+!-----------------------------------------------------------------------
+!
+! Interpolate
+!
+!-----------------------------------------------------------------------
+
+ ! First, compute second derivative values at each node, y2ndDer.
+ call mpas_cubic_spline_coefficients(x,y,n,y2ndDer)
+
+ ! Compute interpolated values yOut.
+ call mpas_interpolate_cubic_spline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,1)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = [',y,'];'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+
+ ! Compute interpolated values yOut.
+ call mpas_integrate_column_cubic_spline( &
+ x,y,y2ndDer,n, &
+ xOut,yOut,nOut)
+
+ ! The following output can be copied directly into Matlab
+ print *, 'subplot(2,1,2)'
+ print '(a,10f8.4,a)', 'x = [',x,'];'
+ print '(a,10f8.4,a)', 'y = 2*cos(-3/2) -2*cos(x/2);'
+ print '(a,100f8.4,a)', 'xOut = [',xOut,'];'
+ print '(a,100f8.4,a)', 'yOut = [',yOut,'];'
+ print *, "plot(x,y,'-*r',xOut,yOut,'x')"
+
+ end subroutine mpas_test_interpolate
+
+end module spline_interpolation
+
Copied: branches/source_renaming/src/operators/mpas_vector_reconstruction.F (from rev 1087, branches/source_renaming/src/operators/module_vector_reconstruction.F)
===================================================================
--- branches/source_renaming/src/operators/mpas_vector_reconstruction.F (rev 0)
+++ branches/source_renaming/src/operators/mpas_vector_reconstruction.F 2011-10-17 17:18:33 UTC (rev 1089)
@@ -0,0 +1,196 @@
+module vector_reconstruction
+
+ use grid_types
+ use configure
+ use constants
+ use rbf_interpolation
+
+ implicit none
+
+ public :: mpas_init_reconstruct, mpas_reconstruct
+
+ contains
+
+ subroutine mpas_init_reconstruct(grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: pre-compute coefficients used by the reconstruct() routine
+ !
+ ! Input: grid meta data
+ !
+ ! Output: grid % coeffs_reconstruct - coefficients used to reconstruct
+ ! velocity vectors at cell centers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: grid
+
+ ! temporary arrays needed in the (to be constructed) init procedure
+ integer :: nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: i, iCell, iEdge, pointCount, maxEdgeCount
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
+ real (kind=RKIND) :: r, cellCenter(3), alpha, tangentPlane(2,3)
+ real (kind=RKIND), allocatable, dimension(:,:) :: edgeOnCellLocations, edgeOnCellNormals, &
+ coeffs
+
+ real(kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors
+ real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ !========================================================
+ ! arrays filled and saved during init procedure
+ !========================================================
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ !========================================================
+ ! temporary variables needed for init procedure
+ !========================================================
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ xEdge => grid % xEdge % array
+ yEdge => grid % yEdge % array
+ zEdge => grid % zEdge % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+ edgeNormalVectors => grid % edgeNormalVectors % array
+ cellTangentPlane => grid % cellTangentPlane % array
+
+
+ ! init arrays
+ coeffs_reconstruct = 0.0
+
+ maxEdgeCount = maxval(nEdgesOnCell)
+
+ allocate(edgeOnCellLocations(maxEdgeCount,3))
+ allocate(edgeOnCellNormals(maxEdgeCount,3))
+ allocate(coeffs(maxEdgeCount,3))
+
+ ! loop over all cells to be solved on this block
+ do iCell=1,nCellsSolve
+ pointCount = nEdgesOnCell(iCell)
+ cellCenter(1) = xCell(iCell)
+ cellCenter(2) = yCell(iCell)
+ cellCenter(3) = zCell(iCell)
+
+ do i=1,pointCount
+ iEdge = edgesOnCell(i,iCell)
+ edgeOnCellLocations(i,1) = xEdge(iEdge)
+ edgeOnCellLocations(i,2) = yEdge(iEdge)
+ edgeOnCellLocations(i,3) = zEdge(iEdge)
+ edgeOnCellNormals(i,:) = edgeNormalVectors(:, iEdge)
+ end do
+
+ alpha = 0.0
+ do i=1,pointCount
+ r = sqrt(sum((cellCenter - edgeOnCellLocations(i,:))**2))
+ alpha = alpha + r
+ enddo
+ alpha = alpha/pointCount
+
+ tangentPlane(1,:) = cellTangentPlane(:,1,iCell)
+ tangentPlane(2,:) = cellTangentPlane(:,2,iCell)
+
+ call mpas_rbf_interp_func_3d_plane_vec_const_dir_comp_coeffs(pointCount, &
+ edgeOnCellLocations(1:pointCount,:), &
+ edgeOnCellNormals(1:pointCount,:), &
+ cellCenter, alpha, tangentPlane, coeffs(1:pointCount,:))
+
+ do i=1,pointCount
+ coeffs_reconstruct(:,i,iCell) = coeffs(i,:)
+ end do
+
+ enddo ! iCell
+
+ deallocate(edgeOnCellLocations)
+ deallocate(edgeOnCellNormals)
+ deallocate(coeffs)
+
+ end subroutine mpas_init_reconstruct
+
+ subroutine mpas_reconstruct(grid, u, uReconstructX, uReconstructY, uReconstructZ, uReconstructZonal, uReconstructMeridional)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Purpose: reconstruct vector field at cell centers based on radial basis functions
+ !
+ ! Input: grid meta data and vector component data residing at cell edges
+ !
+ ! Output: reconstructed vector field (measured in X,Y,Z) located at cell centers
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), dimension(:,:), intent(in) :: u
+ real (kind=RKIND), dimension(:,:), intent(out) :: uReconstructX, uReconstructY, uReconstructZ
+ real (kind=RKIND), dimension(:,:), intent(out) :: uReconstructZonal, uReconstructMeridional
+
+ ! temporary arrays needed in the compute procedure
+ integer :: nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: iCell,iEdge, i
+ real(kind=RKIND), dimension(:), pointer :: latCell, lonCell
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ logical :: on_a_sphere
+
+ real (kind=RKIND) :: clat, slat, clon, slon
+
+
+ ! stored arrays used during compute procedure
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
+
+ ! temporary variables
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+
+ latCell => grid % latCell % array
+ lonCell => grid % lonCell % array
+ on_a_sphere = grid % on_a_sphere
+
+ ! init the intent(out)
+ uReconstructX = 0.0
+ uReconstructY = 0.0
+ uReconstructZ = 0.0
+
+ ! loop over cell centers
+ do iCell=1,nCellsSolve
+ ! a more efficient reconstruction where rbf_values*matrix_reconstruct has been precomputed
+ ! in coeffs_reconstruct
+ do i=1,nEdgesOnCell(iCell)
+ iEdge = edgesOnCell(i,iCell)
+ uReconstructX(:,iCell) = uReconstructX(:,iCell) &
+ + coeffs_reconstruct(1,i,iCell) * u(:,iEdge)
+ uReconstructY(:,iCell) = uReconstructY(:,iCell) &
+ + coeffs_reconstruct(2,i,iCell) * u(:,iEdge)
+ uReconstructZ(:,iCell) = uReconstructZ(:,iCell) &
+ + coeffs_reconstruct(3,i,iCell) * u(:,iEdge)
+
+ enddo
+ enddo ! iCell
+
+ if(on_a_sphere) then
+ do iCell=1,nCellsSolve
+ clat = cos(latCell(iCell))
+ slat = sin(latCell(iCell))
+ clon = cos(lonCell(iCell))
+ slon = sin(lonCell(iCell))
+ uReconstructZonal(:,iCell) = -uReconstructX(:,iCell)*slon + uReconstructY(:,iCell)*clon
+ uReconstructMeridional(:,iCell) = -(uReconstructX(:,iCell)*clon &
+ + uReconstructY(:,iCell)*slon)*slat &
+ + uReconstructZ(:,iCell)*clat
+ end do
+ else
+ uReconstructZonal = uReconstructX
+ uReconstructMeridional = uReconstructY
+ end if
+
+ end subroutine mpas_reconstruct
+
+end module vector_reconstruction
</font>
</pre>