<p><b>xylar@lanl.gov</b> 2010-05-19 11:59:05 -0600 (Wed, 19 May 2010)</p><p>BRANCH COMMIT<br>
<br>
Continued work on some of the Matlab test files<br>
<br>
Removed references to the tanNeu type reconstruction: this appears not to be a numerically feasible way to reconstruct a vector function, at least in my tests so far.<br>
<br>
Removed use of the interp_tests module from core_ocean/module_mpas_interface.F<br>
<br>
Removed redundant variable edgeLocations from the Registry: this variable is already available in the edge{X,Y,Z} variable<br>
<br>
Rewrote the code that computes edge normals and vertical normal vectors to handle lateral boundaries.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/matlab/evaluateFunctionAtDestination.m 2010-05-19 17:59:05 UTC (rev 281)
@@ -1,71 +1,12 @@
-function functionValue = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, functionSamples)
+function [functionValue] = evaluateFunctionAtDestination(matrix, rhs, functionSamples, dim)
-% pointCount = size(sourcePoints,1);
-% matrixSize = pointCount+3;
-% matrix = zeros(matrixSize, matrixSize);
-% rhs = zeros(matrixSize,3);
-%
-% destPoint = sourcePoints(15,:);
-% alpha = 1.35;
-% for(j=1:pointCount)
-% for(i=1:pointCount)
-% rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
-% dotProduct = sum(normals(i,:).*normals(j,:));
-% matrix(i,j) = rbf*dotProduct;
-% end
-% rbf = 1/sqrt(1 + sum((destPoint-sourcePoints(j,:)).^2)/alpha^2);
-% rhs(j,:) = rbf*normals(j,:);
-% matrix(j,pointCount+(1:3)) = normals(j,:);
-% matrix(pointCount+(1:3),j) = normals(j,:);
-% end
-%
-% for(i=1:3)
-% rhs(pointCount+i,i) = 1;
-% end
-%
-% coeffs = matrix\rhs;
-%
-% functionValue(1) = sum(coeffs(1:pointCount,1).*noSlipSamples);
-% functionValue(2) = sum(coeffs(1:pointCount,2).*noSlipSamples);
-% functionValue(3) = sum(coeffs(1:pointCount,3).*noSlipSamples);
+pointCount = length(functionSamples);
-pointCount = size(sourcePoints,1);
-matrixSize = pointCount+3;
-matrix = zeros(matrixSize, matrixSize);
-rhs = zeros(matrixSize,3);
-functionValue = zeros(1,3);
+coeffs = matrix\rhs;
-alpha = 1.35;
-for(j=1:pointCount)
- if(isTangent(j))
- for(i=1:pointCount)
- rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
- normalVector = normals(normalIndices(j),:);
- normalDotX = sum(normalVector.*(sourcePoints(j,:)-sourcePoints(i,:)));
- rbfDerivOverR = -rbf^3;
- dotProduct = sum(normals(i,:).*normals(j,:));
- matrix(i,j) = rbfDerivOverR*dotProduct*normalDotX;
- end
- else
- for(i=1:pointCount)
- rbf = 1/sqrt(1 + sum((sourcePoints(i,:)-sourcePoints(j,:)).^2)/alpha^2);
- dotProduct = sum(normals(i,:).*normals(j,:));
- matrix(i,j) = rbf*dotProduct;
- end
- matrix(j,pointCount+(1:3)) = normals(j,:);
- matrix(pointCount+(1:3),j) = normals(j,:);
- end
- rbf = 1/sqrt(1 + sum((destPoint-sourcePoints(j,:)).^2)/alpha^2);
- rhs(j,:) = rbf*normals(j,:);
-end
+functionValue = zeros(dim,1);
-for(i=1:3)
- rhs(pointCount+i,i) = 1;
+for(i=1:dim)
+ functionValue(i) = sum(coeffs(1:pointCount,i).*functionSamples);
end
-
-coeffs = matrix\rhs;
-
-functionValue(1) = sum(coeffs(1:pointCount,1).*functionSamples);
-functionValue(2) = sum(coeffs(1:pointCount,2).*functionSamples);
-functionValue(3) = sum(coeffs(1:pointCount,3).*functionSamples);
\ No newline at end of file
Modified: branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test_simple.m 2010-05-19 17:59:05 UTC (rev 281)
@@ -2,10 +2,10 @@
close all;
sourcePoints = [-0.5, 0, 0; 0, -.5, 0; 0, .5, 0; .5, 0, 0; .5, 0, 0; .5, 0, 0; 0, 0, -.5; 0, 0, .5];
-%normals = [-1, 0, 0; 0, -1, 0; 0, 1, 0; 1, 0, 0; 0, 1, 0; 0, 0, 1; 0, 0, -1; 0, 0, 1];
-sourcePoints = sourcePoints([1:4,7:8],:);
-normals = zeros(size(sourcePoints));
-normals(:,2) = 1;
+normals = [-1, 0, 0; 0, -1, 0; 0, 1, 0; 1, 0, 0; 0, 1, 0; 0, 0, 1; 0, 0, -1; 0, 0, 1];
+%sourcePoints = sourcePoints([1:4,7:8],:);
+%normals = zeros(size(sourcePoints));
+%normals(:,2) = 1;
normalIndices = zeros(size(sourcePoints,1),1);
isTangent = false(size(normalIndices));
@@ -14,7 +14,7 @@
isTangent(5) = true;
isTangent(6) = true;
-%quiver3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), normals(:,1), normals(:,2), normals(:,3)); axis equal;
+quiver3(sourcePoints(:,1), sourcePoints(:,2), sourcePoints(:,3), normals(:,1), normals(:,2), normals(:,3)); axis equal;
freeSlipFunction = zeros(size(sourcePoints));
freeSlipFunction(:,2) = 1 + sourcePoints(:,1).^2;
@@ -28,10 +28,17 @@
functionValue = zeros(size(X,1),size(X,2),3);
+matrix = constructRBFMatrix(sourcePoints, normals, normalIndices, isTangent, 3);
+
+null(matrix)
+
+return
+
for(yIndex = 1:size(X,2))
for(xIndex = 1:size(X,1))
destPoint = [X(xIndex,yIndex),Y(xIndex,yIndex),0];
- functionValue(xIndex,yIndex,:) = evaluateFunctionAtDestination(destPoint, sourcePoints, normals, normalIndices, isTangent, freeSlipSamples);
+ rhs = constructRBF_RHS(destPoint, sourcePoints, normals, 3);
+ functionValue(xIndex,yIndex,:) = evaluateFunctionAtDestination(matrix, rhs, freeSlipSamples, 3);
end
end
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/Registry 2010-05-19 17:59:05 UTC (rev 281)
@@ -71,8 +71,7 @@
var real areaTriangle ( nVertices ) iro areaTriangle - -
var real edgeNormalVectors ( R3 nEdges ) o edgeNormalVectors - -
-var real radialUnitVectors ( R3 nCells ) o radialUnitVectors - -
-var real edgeLocations ( R3 nEdges ) o edgeLocations - -
+var real localVerticalUnitVectors ( R3 nCells ) o localVerticalUnitVectors - -
var real cellTangentPlane ( R3 TWO nEdges ) o cellTangentPlane - -
var integer cellsOnCell ( maxEdges nCells ) iro cellsOnCell - -
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-05-19 17:59:05 UTC (rev 281)
@@ -18,7 +18,6 @@
use time_integration
use RBF_interpolation
use vector_reconstruction
- use interp_tests
implicit none
@@ -30,12 +29,8 @@
call rbfInterp_initialize(mesh)
- !call init_reconstruct(mesh)
+ call init_reconstruct(mesh)
- call new_init_reconstruct(mesh)
-
- call doTests(mesh)
-
call reconstruct(block % time_levs(1) % state, mesh)
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/Makefile
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/Makefile 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/Makefile 2010-05-19 17:59:05 UTC (rev 281)
@@ -1,6 +1,6 @@
.SUFFIXES: .F .o
-OBJS = module_RBF_interpolation.o module_interp_tests.o module_vector_reconstruction.o
+OBJS = module_RBF_interpolation.o module_vector_reconstruction.o
all: operators
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_RBF_interpolation.F 2010-05-19 17:59:05 UTC (rev 281)
@@ -87,9 +87,9 @@
! 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
+ rbfInterp_func_3DPlane_vec_const_dir_compCoeffs!, &
+ !rbfInterp_func_3D_vec_const_tanNeu_compCoeffs, &
+ !rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
contains
@@ -117,9 +117,8 @@
integer :: nCells, nEdges
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell
integer :: iEdge, iCell, cell1, cell2
- real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
- real(kind=RKIND), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &
- edgeLocations
+ 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
@@ -127,46 +126,53 @@
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
- radialUnitVectors => grid % radialUnitVectors % array
+ localVerticalUnitVectors => grid % localVerticalUnitVectors % array
edgeNormalVectors => grid % edgeNormalVectors % array
- edgeLocations => grid % edgeLocations % array
cellTangentPlane => grid % cellTangentPlane % array
! init arrays
edgeNormalVectors = 0
- radialUnitVectors = 0
+ localVerticalUnitVectors = 0
! loop over all cells to be solved on this block
do iCell=1,nCells
- radialUnitVectors(1,iCell) = xCell(iCell)
- radialUnitVectors(2,iCell) = yCell(iCell)
- radialUnitVectors(3,iCell) = zCell(iCell)
- call unit_vec_in_R3(radialUnitVectors(:,iCell))
+ ! 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
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- ! the normal points from cell1 to cell2
- edgeNormalVectors(1,iEdge) = xCell(cell2) - xCell(cell1)
- edgeNormalVectors(2,iEdge) = yCell(cell2) - yCell(cell1)
- edgeNormalVectors(3,iEdge) = zCell(cell2) - zCell(cell1)
+ 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))
- edgeLocations(1,iEdge) = 0.5*(xCell(cell2) + xCell(cell1))
- edgeLocations(2,iEdge) = 0.5*(yCell(cell2) + yCell(cell1))
- edgeLocations(3,iEdge) = 0.5*(zCell(cell2) + zCell(cell1))
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 = radialUnitVectors(:,iCell)
+ rHat = localVerticalUnitVectors(:,iCell)
normalDotRHat = sum(edgeNormalVectors(:,iEdge)*rHat)
xHatPlane = edgeNormalVectors(:,iEdge) - normalDotRHat*rHat
call unit_vec_in_R3(xHatPlane)
@@ -1185,224 +1191,224 @@
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)
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! ! Purpose: Compute interpolation coefficients in 3D that can be used to
+! ! interpolate a number of vector functions at a given locations. This is useful
+! ! if the interpolation location does not change with time, or if several
+! ! fields are to be interpolated at a given time step. (If both the vector fields
+! ! and the interpolation locations vary with time, there is no clear advantage in
+! ! using either this method or the method for 2D interpoaltion above; for simplicity
+! ! and because we foresee more uses for the method of this subroutine, we have not
+! ! implemented a 3D version of the fixed field, variable interpolation location method
+! ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
+! ! Neumann tangential boundary conditions (such as free slip). The interpolation is
+! ! performed with basis functions that are RBFs plus a constant.
+! ! Input:
+! ! pointCount - the number of "source" points and functionValues supplied
+! ! sourcePoints - the location of the "source" points in the 3D space where the values of
+! ! the function are known
+! ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
+! ! tangent to the interface where the boundary condition will be applied. A Neumann
+! ! boundary condition will be applied at these points in these directions.
+! ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
+! ! gives the normal vector at the same sourcePoint. This information is needed to compute
+! ! the Neumann boundary condition at this point.
+! ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
+! ! is performed by supplying the value of the vector function dotted into each of these unit
+! ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
+! ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
+! ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
+! ! condition and the Neumann tangential boundary conditions at these points.
+! ! destinationPoint - the point where the interpolation will be performed
+! ! alpha - a constant that give the characteristic length scale of the RBFs,
+! ! should be on the order of the distance between points
+! ! Output:
+! ! coefficients - the coefficients used to interpolate a function with Dirichlet
+! ! boundary conditions to the specified destinationPoint
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs(pointCount, &
+! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+! alpha, coefficients)
+!
+! integer, intent(in) :: pointCount
+! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+! logical, dimension(pointCount), intent(in) :: isTangentToInterface
+! integer, dimension(pointCount), intent(in) :: normalVectorIndex
+! real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+! real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+! real(kind=RKIND), intent(in) :: alpha
+! real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!
+! integer :: i, j
+! integer :: matrixSize
+!
+! real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
+! real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+! integer, dimension(:), pointer :: pivotIndices
+!
+! matrixSize = pointCount+3 ! extra space for constant vector
+!
+! allocate(matrix(matrixSize,matrixSize))
+! allocate(matrixCopy(matrixSize,matrixSize))
+! allocate(rhs(matrixSize,3))
+! allocate(coeffs(matrixSize,3))
+! allocate(pivotIndices(matrixSize))
+!
+! matrix = 0.0
+! rhs = 0.0
+! coeffs = 0.0
+!
+! call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
+! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
+! alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+!
+! do i = 1, pointCount
+! if(.not. isTangentToInterface(i)) then
+! matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
+! end if
+! matrix(pointCount+1:pointCount+3,i) &
+! = matrix(i,pointCount+1:pointCount+3)
+! end do
+! do i = 1, 3
+! rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+! end do
+!
+! ! solve each linear system
+! do i = 1, 3
+! matrixCopy = matrix
+! call LEGS(matrixCopy, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+! end do
+! coefficients = coeffs(1:pointCount,:)
+!
+! deallocate(matrix)
+! deallocate(matrixCopy)
+! deallocate(rhs)
+! deallocate(coeffs)
+! deallocate(pivotIndices)
+!
+! end subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs
+!
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! ! Purpose: Compute interpolation coefficients in 3D that can be used to
+! ! interpolate a number of vector functions at a given locations. This is useful
+! ! if the interpolation location does not change with time, or if several
+! ! fields are to be interpolated at a given time step. (If both the vector fields
+! ! and the interpolation locations vary with time, there is no clear advantage in
+! ! using either this method or the method for 2D interpoaltion above; for simplicity
+! ! and because we foresee more uses for the method of this subroutine, we have not
+! ! implemented a 3D version of the fixed field, variable interpolation location method
+! ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
+! ! Neumann tangential boundary conditions (such as free slip). The interpolation is
+! ! performed with basis functions that are RBFs plus a constant.
+! ! Input:
+! ! pointCount - the number of "source" points and functionValues supplied
+! ! sourcePoints - the location of the "source" points in the 3D space where the values of
+! ! the function are known. The sourcePoints are projected into the plane given by
+! ! planeBasisVectors
+! ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
+! ! tangent to the interface where the boundary condition will be applied. A Neumann
+! ! boundary condition will be applied at these points in these directions.
+! ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
+! ! gives the normal vector at the same sourcePoint. This information is needed to compute
+! ! the Neumann boundary condition at this point.
+! ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
+! ! is performed by supplying the value of the vector function dotted into each of these unit
+! ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
+! ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
+! ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
+! ! condition and the Neumann tangential boundary conditions at these points. The unitVectors
+! ! are projected into the plane given by planeBasisVectors
+! ! destinationPoint - the point where the interpolation will be performed. The destinationPoint
+! ! is projected into the plane given by planeBasisVectors
+! ! alpha - a constant that give the characteristic length scale of the RBFs,
+! ! should be on the order of the distance between points
+! ! planeBasisVectors - the basis fectors for the plane where interpolation is performed.
+! ! All points are projected into this plane.
+! ! Output:
+! ! coefficients - the coefficients used to interpolate a function with Dirichlet
+! ! boundary conditions to the specified destinationPoint
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs(&
+! pointCount, sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, &
+! destinationPoint, alpha, planeBasisVectors, coefficients)
+!
+! integer, intent(in) :: pointCount
+! real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
+! logical, dimension(pointCount), intent(in) :: isTangentToInterface
+! integer, dimension(pointCount), intent(in) :: normalVectorIndex
+! real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
+! real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
+! real(kind=RKIND), intent(in) :: alpha
+! real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+! real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
+!
+! integer :: i, j
+! integer :: matrixSize
+!
+! real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
+! real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
+! real(kind=RKIND), dimension(2) :: planarDestinationPoint
+!
+! real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
+! real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
+! integer, dimension(:), pointer :: pivotIndices
+!
+! matrixSize = pointCount+2 ! space for constant vector in plane
+!
+! allocate(matrix(matrixSize,matrixSize))
+! allocate(matrixCopy(matrixSize,matrixSize))
+! allocate(rhs(matrixSize,2))
+! allocate(coeffs(matrixSize,2))
+! allocate(pivotIndices(matrixSize))
+!
+! matrix = 0.0
+! rhs = 0.0
+! coeffs = 0.0
+!
+! do i = 1, pointCount
+! planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
+! planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
+! planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+! planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+! end do
+! planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
+! planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
+! call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
+! planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
+! planarDestinationPoint, alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
+!
+! do i = 1, pointCount
+! if(.not. isTangentToInterface(i)) then
+! matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
+! matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+! end if
+! matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
+! end do
+! do i = 1,2
+! rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
+! end do
+!
+! ! solve each linear system
+! matrixCopy = matrix
+! call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
+! call LEGS(matrixCopy, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
+!
+! coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
+! + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
+! coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
+! + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
+! coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
+! + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
+! coefficients = coeffs(1:pointCount,:)
+!
+! deallocate(matrix)
+! deallocate(matrixCopy)
+! deallocate(rhs)
+! deallocate(coeffs)
+! deallocate(pivotIndices)
+!
+! end subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isTangentToInterface
- integer, dimension(pointCount), intent(in) :: normalVectorIndex
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
- real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
- integer, dimension(:), pointer :: pivotIndices
-
- matrixSize = pointCount+3 ! extra space for constant vector
-
- allocate(matrix(matrixSize,matrixSize))
- allocate(matrixCopy(matrixSize,matrixSize))
- allocate(rhs(matrixSize,3))
- allocate(coeffs(matrixSize,3))
- allocate(pivotIndices(matrixSize))
-
- matrix = 0.0
- rhs = 0.0
- coeffs = 0.0
-
- call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
- sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
- alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-
- do i = 1, pointCount
- if(.not. isTangentToInterface(i)) then
- matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
- end if
- matrix(pointCount+1:pointCount+3,i) &
- = matrix(i,pointCount+1:pointCount+3)
- end do
- do i = 1, 3
- rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
- end do
-
- ! solve each linear system
- do i = 1, 3
- matrixCopy = matrix
- call LEGS(matrixCopy, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
- end do
- coefficients = coeffs(1:pointCount,:)
-
- deallocate(matrix)
- deallocate(matrixCopy)
- deallocate(rhs)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- end subroutine rbfInterp_func_3D_vec_const_tanNeu_compCoeffs
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Purpose: Compute interpolation coefficients in 3D that can be used to
- ! interpolate a number of vector functions at a given locations. This is useful
- ! if the interpolation location does not change with time, or if several
- ! fields are to be interpolated at a given time step. (If both the vector fields
- ! and the interpolation locations vary with time, there is no clear advantage in
- ! using either this method or the method for 2D interpoaltion above; for simplicity
- ! and because we foresee more uses for the method of this subroutine, we have not
- ! implemented a 3D version of the fixed field, variable interpolation location method
- ! as we have in 2D.) Coefficients are produced for handling Dirichlet normal /
- ! Neumann tangential boundary conditions (such as free slip). The interpolation is
- ! performed with basis functions that are RBFs plus a constant.
- ! Input:
- ! pointCount - the number of "source" points and functionValues supplied
- ! sourcePoints - the location of the "source" points in the 3D space where the values of
- ! the function are known. The sourcePoints are projected into the plane given by
- ! planeBasisVectors
- ! isTangentToInterface - a logical array indicating which sourcePoints/unitVectors are
- ! tangent to the interface where the boundary condition will be applied. A Neumann
- ! boundary condition will be applied at these points in these directions.
- ! normalVectorIndex - where isTangentToInterface == .true., the index into unitVectors that
- ! gives the normal vector at the same sourcePoint. This information is needed to compute
- ! the Neumann boundary condition at this point.
- ! unitVectors - the unit vectors associated with each of the sourcePoints. Interpolation
- ! is performed by supplying the value of the vector function dotted into each of these unit
- ! vectors. If multiple unit vectors are supplied at the same sourcePoint, they *must* be
- ! orthogonal for the interpolation to succeed. A normal vector and two tangential vectors
- ! are needed at each interface point in order to satisfy the Dirichlet normal boundary
- ! condition and the Neumann tangential boundary conditions at these points. The unitVectors
- ! are projected into the plane given by planeBasisVectors
- ! destinationPoint - the point where the interpolation will be performed. The destinationPoint
- ! is projected into the plane given by planeBasisVectors
- ! alpha - a constant that give the characteristic length scale of the RBFs,
- ! should be on the order of the distance between points
- ! planeBasisVectors - the basis fectors for the plane where interpolation is performed.
- ! All points are projected into this plane.
- ! Output:
- ! coefficients - the coefficients used to interpolate a function with Dirichlet
- ! boundary conditions to the specified destinationPoint
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs(&
- pointCount, sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, &
- destinationPoint, alpha, planeBasisVectors, coefficients)
-
- integer, intent(in) :: pointCount
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: sourcePoints
- logical, dimension(pointCount), intent(in) :: isTangentToInterface
- integer, dimension(pointCount), intent(in) :: normalVectorIndex
- real(kind=RKIND), dimension(pointCount,3), intent(in) :: unitVectors
- real(kind=RKIND), dimension(3), intent(in) :: destinationPoint
- real(kind=RKIND), intent(in) :: alpha
- real(kind=RKIND), dimension(2,3) :: planeBasisVectors
- real(kind=RKIND), dimension(pointCount, 3), intent(out) :: coefficients
-
- integer :: i, j
- integer :: matrixSize
-
- real(kind=RKIND), dimension(pointCount,2) :: planarSourcePoints
- real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
- real(kind=RKIND), dimension(2) :: planarDestinationPoint
-
- real(kind=RKIND), dimension(:,:), pointer :: matrix, matrixCopy
- real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
- integer, dimension(:), pointer :: pivotIndices
-
- matrixSize = pointCount+2 ! space for constant vector in plane
-
- allocate(matrix(matrixSize,matrixSize))
- allocate(matrixCopy(matrixSize,matrixSize))
- allocate(rhs(matrixSize,2))
- allocate(coeffs(matrixSize,2))
- allocate(pivotIndices(matrixSize))
-
- matrix = 0.0
- rhs = 0.0
- coeffs = 0.0
-
- do i = 1, pointCount
- planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
- planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
- planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
- planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
- end do
- planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
- planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
- call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
- planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
- planarDestinationPoint, alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
-
- do i = 1, pointCount
- if(.not. isTangentToInterface(i)) then
- matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
- matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
- end if
- matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
- end do
- do i = 1,2
- rhs(pointCount+i,i) = 1.0 ! the unit vector in the ith direction
- end do
-
- ! solve each linear system
- matrixCopy = matrix
- call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
- call LEGS(matrixCopy, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
-
- coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
- + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
- coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
- + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
- coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
- + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
- coefficients = coeffs(1:pointCount,:)
-
- deallocate(matrix)
- deallocate(matrixCopy)
- deallocate(rhs)
- deallocate(coeffs)
- deallocate(pivotIndices)
-
- end subroutine rbfInterp_func_3DPlane_vec_const_tanNeu_compCoeffs
-
!!!!!!!!!!!!!!!!!!!!!
! private subroutines
!!!!!!!!!!!!!!!!!!!!!
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-05-19 17:59:05 UTC (rev 281)
@@ -54,29 +54,30 @@
integer :: iCell, iEdge
- real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
- real(kind=RKIND), dimension(:,:), pointer :: radialUnitVectors, edgeNormalVectors, &
- edgeLocations
+ real(kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell, xEdge, yEdge, zEdge
+ real(kind=RKIND), dimension(:,:), pointer :: localVerticalUnitVectors, edgeNormalVectors
real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
- radialUnitVectors => grid % radialUnitVectors % array
+ localVerticalUnitVectors => grid % localVerticalUnitVectors % array
edgeNormalVectors => grid % edgeNormalVectors % array
- edgeLocations => grid % edgeLocations % array
cellTangentPlane => grid % cellTangentPlane % array
xCell => grid % xCell % array
yCell => grid % yCell % array
zCell => grid % zCell % array
+ xEdge => grid % xEdge % array
+ yEdge => grid % yEdge % array
+ zEdge => grid % zEdge % array
iCell = 1
iEdge = grid%edgesOnCell%array(1,iCell)
- print *, "cell: ", iCell, "radial unit vector: ", radialUnitVectors(:,iCell)
+ print *, "cell: ", iCell, "radial unit vector: ", localVerticalUnitVectors(:,iCell)
print *, "x, y, z: ", xCell(iCell), yCell(iCell), zCell(iCell)
print *, "tangentPlane: ", cellTangentPlane(:,:,iCell)
print *, "edge: ", iEdge, "normal: ", edgeNormalVectors(:,iEdge)
- print *, "location: ", edgeLocations(:,iEdge)
+ print *, "location: ", xEdge(iEdge), yEdge(iEdge), zEdge(iEdge)
end subroutine test_ibInterp_initialize
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-05-19 17:25:39 UTC (rev 280)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-05-19 17:59:05 UTC (rev 281)
@@ -1,204 +1,17 @@
module vector_reconstruction
- use grid_types
- use configure
- use constants
- use RBF_interpolation
+ use grid_types
+ use configure
+ use constants
+ use RBF_interpolation
- implicit none
+ implicit none
- public :: init_reconstruct, reconstruct, new_init_reconstruct
+ public :: init_reconstruct, reconstruct
- contains
+ 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 (grid_meta), intent(inout) :: grid
-
- ! temporary arrays needed in the (to be constructed) init procedure
- integer :: nCells, nEdges, nVertLevels, nCellsSolve
- integer, dimension(:,:), pointer :: cellsOnCell, edgesOnCell, cellsOnEdge
- integer, dimension(:), pointer :: nEdgesOnCell
- integer :: iEdge, iCell, k, cell1, cell2, EdgeMax, j, i, npoints, matrixSize
- integer :: lwork, info
- integer, allocatable, dimension(:) :: pivotingIndices
- real (kind=RKIND), dimension(:), pointer :: dcEdge, xCell, yCell, zCell
- real (kind=RKIND) :: r, rbfValue, v, X1(3), X2(3), alpha, rHat(3), normalDotRHat
- real (kind=RKIND) :: xPlane, yPlane, xNormalPlane, yNormalPlane, xHatPlane(3), yHatPlane(3)
- real (kind=RKIND), allocatable, dimension(:,:,:) :: xLoc
-
- real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
- real (kind=RKIND), allocatable, dimension(:,:) :: mwork
- real (kind=RKIND), dimension(:,:), pointer :: matrix, invMatrix
- real (kind=RKIND), dimension(:,:), pointer :: normals
-
- !========================================================
- ! 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
- cellsOnCell => grid % cellsOnCell % array
- cellsOnEdge => grid % cellsOnEdge % array
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- dcEdge => grid % dcEdge % array
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
-
- ! allocate arrays
- EdgeMax = maxval(nEdgesOnCell)
- allocate(xLoc(3,EdgeMax,nCells))
-
- allocate(normals(3,EdgeMax))
-
- ! init arrays
- coeffs_reconstruct = 0.0
- normals = 0
-
- ! loop over all cells to be solved on this block
- do iCell=1,nCellsSolve
-
- ! fill normal vector and xLoc arrays
- ! X1 is the location of the reconstruction, X2 are the neighboring centers,
- ! xLoc is the edge positions
- cell1 = iCell
- X1(1) = xCell(cell1)
- X1(2) = yCell(cell1)
- X1(3) = zCell(cell1)
-
- rHat = X1
- call unit_vector_in_R3(rHat)
-
- do j=1,nEdgesOnCell(iCell)
- iEdge = edgesOnCell(j,iCell)
- if (iCell == cellsOnEdge(1,iEdge)) then
- cell2 = cellsOnEdge(2,iEdge)
- X2(1) = xCell(cell2)
- X2(2) = yCell(cell2)
- X2(3) = zCell(cell2)
- normals(:,j) = X2(:) - X1(:)
- else
- cell2 = cellsOnEdge(1,iEdge)
- X2(1) = xCell(cell2)
- X2(2) = yCell(cell2)
- X2(3) = zCell(cell2)
- normals(:,j) = X1(:) - X2(:)
- endif
-
-
- call unit_vector_in_R3(normals(:,j))
-
- xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
- enddo
-
- npoints = nEdgesOnCell(iCell) ! only loop over the number of edges for this cell
- matrixSize = npoints+2 ! room for 2 vector components for constant flow
- ! basis functions
- allocate(matrix(matrixSize,matrixSize))
- matrix = 0.0
- alpha = 0.0
- do i=1,npoints
- call get_distance(xLoc(:,i,iCell),X1(:),r)
- alpha = alpha + r
- enddo
- alpha = alpha / npoints
- do j=1,npoints
- do i=1,npoints
- call get_distance(xLoc(:,i,iCell),xLoc(:,j,iCell),r)
- r = r / alpha
- call evaluate_rbf(r,rbfValue)
- call get_dotproduct(normals(:,i),normals(:,j),v)
- matrix(i,j) = v*rbfValue
- enddo
- enddo
-
-
- ! add matrix entries to suppoert constant flow
- ! xHat and yHat are a local basis in the plane of the horizontal cell
- ! we arbitrarily choose xHat to point toward the first edge
- call get_dotproduct(normals(:,1),rHat,normalDotRHat)
- xHatPlane = normals(:,1) - normalDotRHat*rHat(:)
- call unit_vector_in_R3(xHatPlane)
-
- call cross_product_in_R3(rHat, xHatPlane, yHatPlane)
- call unit_vector_in_R3(yHatPlane) ! just to be sure...
-
- do j=1,npoints
- call get_dotproduct(normals(:,j),xHatPlane, xNormalPlane)
- call get_dotproduct(normals(:,j),yHatPlane, yNormalPlane)
- matrix(j,npoints+1) = xNormalPlane
- matrix(j,npoints+2) = yNormalPlane
- matrix(npoints+1,j) = matrix(j,npoints+1)
- matrix(npoints+2,j) = matrix(j,npoints+2)
- end do
-
- ! invert matrix
- allocate(invMatrix(matrixSize,matrixSize))
- allocate(pivotingIndices(matrixSize))
- invMatrix = 0.0
- pivotingIndices = 0
- call migs(matrix,matrixSize,invMatrix,pivotingIndices)
-
- ! compute the coefficients for reconstructing uX, uY, uZ at cell centers from
- ! u_i normal to edges
- ! uX = sum_j(coeffs(1,j) * u_j) (similarly for Y,Z)
- ! coeffs(:,j) = sum_i(rbf_values(i) * normal(:,i) * matrix(i,j))
- do i=1,npoints
- ! compute value of RBF when evaluated between reconstruction location and edge locations
- call get_distance(xLoc(:,i,iCell), X1(:), r)
- r = r / alpha
- call evaluate_rbf(r,rbfValue)
-
- ! project the normals onto tangent plane of the cell
- ! normal = normal - (normal dot r_hat) r_hat
- ! rHat, the unit vector pointing from the domain center to the cell center
- call get_dotproduct(normals(:,i),rHat,normalDotRHat)
- normals(:,i) = normals(:,i) - normalDotRHat*rHat(:)
-
- do j=1,npoints
- coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
- + rbfValue * normals(:,i) * invMatrix(i,j)
- enddo
- enddo
-
- ! polynomial parts
- do j=1,npoints
- coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
- + invMatrix(npoints+1,j)*xHatPlane
- coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
- + invMatrix(npoints+2,j)*yHatPlane
- enddo
-
- deallocate(matrix)
- deallocate(invMatrix)
- deallocate(pivotingIndices)
-
- enddo ! iCell
-
- deallocate(xLoc)
- deallocate(normals)
-
- end subroutine init_reconstruct
-
- subroutine new_init_reconstruct(grid)
+ subroutine init_reconstruct(grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose: pre-compute coefficients used by the reconstruct() routine
!
@@ -216,13 +29,13 @@
integer :: nCellsSolve
integer, dimension(:,:), pointer :: edgesOnCell
integer, dimension(:), pointer :: nEdgesOnCell
- integer :: i, iCell, pointCount, maxEdgeCount
- real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+ 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, edgeLocations
+ real(kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors
real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
@@ -238,11 +51,13 @@
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
- edgeLocations => grid % edgeLocations % array
cellTangentPlane => grid % cellTangentPlane % array
@@ -263,8 +78,11 @@
cellCenter(3) = zCell(iCell)
do i=1,pointCount
- edgeOnCellLocations(i,:) = edgeLocations(:, edgesOnCell(i,iCell))
- edgeOnCellNormals(i,:) = edgeNormalVectors(:, edgesOnCell(i,iCell))
+ 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
@@ -292,7 +110,7 @@
deallocate(edgeOnCellNormals)
deallocate(coeffs)
- end subroutine new_init_reconstruct
+ end subroutine init_reconstruct
subroutine reconstruct(state, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -352,202 +170,4 @@
end subroutine reconstruct
-
- subroutine get_distance(x1,x2,xout)
- implicit none
- real(kind=RKIND), intent(in) :: x1(3), x2(3)
- real(kind=RKIND), intent(out) :: xout
- xout = sqrt( (x1(1)-x2(1))**2 + (x1(2)-x2(2))**2 + (x1(3)-x2(3))**2 )
- end subroutine get_distance
-
- subroutine get_dotproduct(x1,x2,xout)
- implicit none
- real(kind=RKIND), intent(in) :: x1(3), x2(3)
- real(kind=RKIND), intent(out) :: xout
- xout = x1(1)*x2(1) + x1(2)*x2(2) + x1(3)*x2(3)
- end subroutine get_dotproduct
-
-
- subroutine unit_vector_in_R3(xin)
- implicit none
- real (kind=RKIND), intent(inout) :: xin(3)
- real (kind=RKIND) :: mag
- mag = sqrt(xin(1)**2+xin(2)**2+xin(3)**2)
- xin(:) = xin(:) / mag
- end subroutine unit_vector_in_R3
-
-
- subroutine evaluate_rbf(xin,xout)
- real(kind=RKIND), intent(in) :: xin
- real(kind=RKIND), intent(out) :: xout
-
- ! Gaussian
- ! xout = exp(-r**2)
-
- ! multiquadrics
- xout = 1.0 / sqrt(1.0**2 + xin**2)
-
- ! other
- ! xout = 1.0 / (1.0**2 + r**2)
-
- end subroutine evaluate_rbf
-
-!======================================================================
-! BEGINNING OF CROSS_PRODUCT_IN_R3
-!======================================================================
- subroutine cross_product_in_R3(p_1,p_2,p_out)
-
-!-----------------------------------------------------------------------
-! PURPOSE: compute p_1 cross p_2 and place in p_out
-!-----------------------------------------------------------------------
-
-!-----------------------------------------------------------------------
-! intent(in)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(in) :: &
- p_1 (3), &
- p_2 (3)
-
-!-----------------------------------------------------------------------
-! intent(out)
-!-----------------------------------------------------------------------
- real (kind=RKIND), intent(out) :: &
- p_out (3)
-
- p_out(1) = p_1(2)*p_2(3)-p_1(3)*p_2(2)
- p_out(2) = p_1(3)*p_2(1)-p_1(1)*p_2(3)
- p_out(3) = p_1(1)*p_2(2)-p_1(2)*p_2(1)
-
- end subroutine cross_product_in_R3
-!======================================================================
-! END OF CROSS_PRODUCT_IN_R3
-!======================================================================
-
-! Updated 10/24/2001.
-!
-!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.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 vector_reconstruction
</font>
</pre>