<p><b>xylar@lanl.gov</b> 2010-04-28 22:02:42 -0600 (Wed, 28 Apr 2010)</p><p>Branch Commit<br>
<br>
Fixed several bugs in ibInterp_computeVectorDirichletPlanarRBFCoefficients, the routine that computes interpolation coefficients for a vector quantity from edges that lie in a plane (e.g., in the local tangent plane of a cell).<br>
<br>
Added new_init_reconstruct routine to the vector_reconstruction module that makes use of ibInterp_computeVectorDirichletPlanarRBFCoefficients and of the geometry information computed in ibInterp_initialize. The new routine produces the same coefficients as the old routine to near machine precision.<br>
<br>
Under the ocean core:<br>
Replaced the call to init_reconstruct with a call to new_init_reconstruct. I will delete the old routine when we're sure the new one works.<br>
<br>
Uncommented the call to reconstruct in time_integration and removed some unneeded comments<br>
<br>
Added a call to reconstruct in the mpas_init routine in the mpas_interface so uReconstruct{X,Y,Z} are computed at initialization as well<br>
<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m
===================================================================
--- branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/matlab/vector_test.m 2010-04-29 04:02:42 UTC (rev 223)
@@ -82,4 +82,33 @@
freeSlipFunction(interfaceMask,3) = planeNormal(1);
noSlipSamples = sum(noSlipFunction.*normals,2);
-freeSlipSamples = sum(freeSlipFunction.*normals,2);
\ No newline at end of file
+freeSlipSamples = sum(freeSlipFunction.*normals,2);
+
+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);
\ No newline at end of file
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/module_time_integration.F 2010-04-29 04:02:42 UTC (rev 223)
@@ -4,9 +4,7 @@
use configure
use constants
use dmpar
- ! xsad 10-02-05:
use vector_reconstruction
- ! xsad 10-02-05 end
contains
@@ -225,10 +223,7 @@
call compute_solve_diagnostics(dt, block % time_levs(2) % state, block % mesh)
- ! xsad 10-02-09:
- ! commenting this out until we incorporate the necessary lapack routines into mpas
- !call reconstruct(block % time_levs(2) % state, block % mesh)
- ! xsad 10-02-09 end
+ call reconstruct(block % time_levs(2) % state, block % mesh)
block => block % next
end do
Modified: branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/src/core_ocean/mpas_interface.F 2010-04-29 04:02:42 UTC (rev 223)
@@ -18,9 +18,7 @@
use time_integration
use IB_interpolation
use vector_reconstruction
- use interp_tests
-
implicit none
type (block_type), intent(inout) :: block
@@ -31,10 +29,12 @@
call ibInterp_initialize(mesh)
- call doTests(mesh)
+ !call init_reconstruct(mesh)
- call init_reconstruct(mesh)
+ call new_init_reconstruct(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
! input arguement into mpas_init. Ask about that later. For now, there will be
! no initial statistics write.
@@ -78,10 +78,7 @@
call timestep(domain, dt)
- ! mrp 100120:
if (mod(itimestep, config_stats_interval) == 0) then
- ! xsad 10-02-08:
- !call write_stats(domain, itimestep, dt)
block_ptr => domain % blocklist
if(associated(block_ptr % next)) then
write(0,*) 'Error: computeGlobalDiagnostics assumes ',&
@@ -93,9 +90,7 @@
block_ptr % time_levs(2) % state, block_ptr % mesh, &
itimestep, dt)
call timer_stop("global diagnostics")
- ! xsad 10-02-08 end
end if
- ! mrp 100120 end
end subroutine mpas_timestep
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_IB_interpolation.F 2010-04-29 04:02:42 UTC (rev 223)
@@ -14,8 +14,8 @@
ibInterp_computeScalarRBFAndConstantCoefficients, &
ibInterp_computeScalarRBFAndPlaneCoefficients, &
ibInterp_computeScalarRBFAndLinearCoefficients, &
- ibInterp_computeVectorNoSlipRBFCoefficients, &
- ibInterp_computeVectorNoSlipPlanarRBFCoefficients, &
+ ibInterp_computeVectorDirichletRBFCoefficients, &
+ ibInterp_computeVectorDirichletPlanarRBFCoefficients, &
ibInterp_computeVectorFreeSlipRBFCoefficients, &
ibInterp_computeVectorFreeSlipPlanarRBFCoefficients
@@ -441,7 +441,7 @@
end subroutine ibInterp_computeScalarRBFAndLinearCoefficients
- subroutine ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &
+ subroutine ibInterp_computeVectorDirichletRBFCoefficients(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients)
@@ -455,13 +455,14 @@
integer :: i, j
integer :: matrixSize
- real(kind=RKIND), dimension(:,:), pointer :: matrix
+ 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))
@@ -470,9 +471,9 @@
rhs = 0.0
coeffs = 0.0
- call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 3, &
+ call setUpVectorDirichletRBFMatrixAndRHS(pointCount, 3, &
sourcePoints, unitVectors, destinationPoint, &
- alpha, matrix, rhs)
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
do i = 1, pointCount
matrix(i,pointCount+1:pointCount+3) = unitVectors(i,:)
@@ -485,18 +486,20 @@
! solve each linear system
do i = 1, 3
- call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ 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 ibInterp_computeVectorNoSlipRBFCoefficients
+ end subroutine ibInterp_computeVectorDirichletRBFCoefficients
- subroutine ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &
+ subroutine ibInterp_computeVectorDirichletPlanarRBFCoefficients(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, planeBasisVectors, coefficients)
@@ -517,13 +520,14 @@
real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
real(kind=RKIND), dimension(2) :: planarDestinationPoint
- real(kind=RKIND), dimension(:,:), pointer :: matrix
+ 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))
@@ -540,35 +544,37 @@
end do
planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
- call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 2, &
+
+ call setUpVectorDirichletRBFMatrixAndRHS(pointCount, 2, &
planarSourcePoints, planarUnitVectors, planarDestinationPoint, &
- alpha, matrix, rhs)
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
do i = 1, pointCount
- matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
- matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
+ matrix(i,pointCount+1:pointCount+2) = planarUnitVectors(i,:)
matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
end do
- rhs(pointCount+1,:) = planeBasisVectors(1,:)
- rhs(pointCount+2,:) = planeBasisVectors(2,:)
+ 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(matrix, matrixSize, rhs(:,2), coeffs(:,2), 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)
+ 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 ibInterp_computeVectorNoSlipPlanarRBFCoefficients
+ end subroutine ibInterp_computeVectorDirichletPlanarRBFCoefficients
subroutine ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &
sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
@@ -586,13 +592,14 @@
integer :: i, j
integer :: matrixSize
- real(kind=RKIND), dimension(:,:), pointer :: matrix
+ 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))
@@ -603,7 +610,7 @@
call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 3, &
sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
- alpha, matrix, rhs)
+ alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
do i = 1, pointCount
if(.not. isTangentToInterface(i)) then
@@ -618,11 +625,13 @@
! solve each linear system
do i = 1, 3
- call LEGS(matrix, matrixSize, rhs(:,i), coeffs(:,i), pivotIndices)
+ 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)
@@ -650,13 +659,14 @@
real(kind=RKIND), dimension(pointCount,2) :: planarUnitVectors
real(kind=RKIND), dimension(2) :: planarDestinationPoint
- real(kind=RKIND), dimension(:,:), pointer :: matrix
+ 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))
@@ -675,7 +685,7 @@
planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
call setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, 2, &
planarSourcePoints, isTangentToInterface, normalVectorIndex, planarUnitVectors, &
- planarDestinationPoint, alpha, matrix, rhs)
+ planarDestinationPoint, alpha, matrix(1:pointCount,1:pointCount), rhs(1:pointCount,:))
do i = 1, pointCount
if(.not. isTangentToInterface(i)) then
@@ -684,12 +694,14 @@
end if
matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
end do
- rhs(pointCount+1,:) = planeBasisVectors(1,:)
- rhs(pointCount+2,:) = planeBasisVectors(2,:)
+ 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(matrix, matrixSize, rhs(:,2), coeffs(:,2), 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)
@@ -700,6 +712,7 @@
coefficients = coeffs(1:pointCount,:)
deallocate(matrix)
+ deallocate(matrixCopy)
deallocate(rhs)
deallocate(coeffs)
deallocate(pivotIndices)
@@ -793,7 +806,7 @@
end subroutine setUpScalarRBFMatrixAndRHS
- subroutine setUpVectorNoSlipRBFMatrixAndRHS(pointCount, dimensions, &
+ subroutine setUpVectorDirichletRBFMatrixAndRHS(pointCount, dimensions, &
sourcePoints, unitVectors, destinationPoint, &
alpha, matrix, rhs)
@@ -824,7 +837,7 @@
rhs(j,:) = evaluateRBF(rSquared)*unitVectors(j,:)
end do
- end subroutine setUpVectorNoSlipRBFMatrixAndRHS
+ end subroutine setUpVectorDirichletRBFMatrixAndRHS
subroutine setUpVectorFreeSlipRBFMatrixAndRHS(pointCount, dimensions, &
sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_interp_tests.F 2010-04-29 04:02:42 UTC (rev 223)
@@ -32,9 +32,12 @@
print *, "ScalarRBFAndPlanar test"
call test_ibInterp_computeScalarRBFAndPlane
-print *, "VectorNoSlipRBFCoefficients test"
- call test_ibInterp_computeVectorNoSlipRBFCoefficients
+print *, "VectorDirichletRBFCoefficients test"
+ call test_ibInterp_computeVectorDirichletRBFCoefficients
+print *, "VectorDirichletPlanarRBFCoefficients test"
+ call test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
+
print *, "tests finished, stopping..."
stop
@@ -345,8 +348,56 @@
print *, "deriv should be:", testFunction(12)
end subroutine test_ibInterp_computeScalarRBFAndPlane
- subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients
+ subroutine test_ibInterp_computeVectorDirichletRBFCoefficients2
+ integer, parameter :: pointCount = 6
+ real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,3) :: unitVectors
+ real(kind=RKIND), dimension(3) :: destinationPoint, functionValue
+ real(kind=RKIND) :: alpha
+ real(kind=RKIND), dimension(pointCount, 3) :: coefficients
+
+ real(kind=RKIND), dimension(pointCount) :: testFunction
+
+ alpha = 1.35
+
+ sourcePoints(1,:) = (/0.0, 0.0, -0.500000000000000/)
+ sourcePoints(2,:) = (/0.0, -0.500000000000000, 0.0/)
+ sourcePoints(3,:) = (/-0.500000000000000, 0.0, 0.0/)
+ sourcePoints(4,:) = (/0.500000000000000, 0.0, 0.0/)
+ sourcePoints(5,:) = (/0.0, 0.500000000000000, 0.0/)
+ sourcePoints(6,:) = (/0.0, 0.0, 0.500000000000000/)
+
+ unitVectors(1,:) = (/0.0, 0.0, 1.0/)
+ unitVectors(2,:) = (/0.0, 1.0, 0.0/)
+ unitVectors(3,:) = (/1.0, 0.0, 0.0/)
+ unitVectors(4,:) = (/1.0, 0.0, 0.0/)
+ unitVectors(5,:) = (/0.0, 1.0, 0.0/)
+ unitVectors(6,:) = (/0.0, 0.0, 1.0/)
+
+ testFunction = (/ -0.500000000000000, -0.500000000000000, -0.500000000000000, 0.500000000000000, 0.500000000000000, 0.500000000000000 /)
+
+ destinationPoint = sourcePoints(3,:)
+
+ call ibInterp_computeVectorDirichletRBFCoefficients(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, coefficients)
+
+print *, "coefficients", coefficients
+
+ functionValue(1) = sum(coefficients(:,1)*testFunction)
+ functionValue(2) = sum(coefficients(:,2)*testFunction)
+ functionValue(3) = sum(coefficients(:,3)*testFunction)
+
+print *, "function value:", functionValue
+print *, "function should be approximately:", sourcePoints(3,1), sourcePoints(3,2), sourcePoints(3,3)
+print *, "function value dot normal:", sum(functionValue*unitVectors(3,:))
+print *, "function dot normal should be:", testFunction(3)
+
+ end subroutine test_ibInterp_computeVectorDirichletRBFCoefficients2
+
+ subroutine test_ibInterp_computeVectorDirichletRBFCoefficients
+
integer, parameter :: pointCount = 46
real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
real(kind=RKIND), dimension(pointCount,3) :: unitVectors
@@ -454,97 +505,94 @@
testFunction = (/-1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 1.0, 0.0, -0.250000000000000, 0.0, -0.250000000000000, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.822383489827826, 0.864313506097250, 0.250000000000000, 0.271157178383451, -0.0724644183610947, 0.448422287815161, -0.0506471741233457, 0.582052908771373, 0.106212968041363, 0.185460903190116, -0.624037720753383, -0.367545123997483, -0.277350098112615, 0.0881719303783468, -0.416025147168922, 0.0352255506168395, 0.138675049056307, -0.125126342375154, 0.762712769809690, -0.0343611199268557, 0.970725343394151, -0.199251274041264, 0.277350098112615, -0.0247308217712550, 0.485362671697076/)
- destinationPoint = sourcePoints(15,:)
+ destinationPoint = sourcePoints(10,:)
- call ibInterp_computeVectorNoSlipRBFCoefficients(pointCount, &
+ call ibInterp_computeVectorDirichletRBFCoefficients(pointCount, &
sourcePoints, unitVectors, destinationPoint, &
alpha, coefficients)
+print *, "coefficients", coefficients
+
functionValue(1) = sum(coefficients(:,1)*testFunction)
- functionValue(2) = sum(coefficients(:,1)*testFunction)
- functionValue(3) = sum(coefficients(:,1)*testFunction)
+ functionValue(2) = sum(coefficients(:,2)*testFunction)
+ functionValue(3) = sum(coefficients(:,3)*testFunction)
print *, "function value:", functionValue
-print *, "function value dot normal:", sum(functionValue*unitVectors(15,:))
-print *, "function should be:", testFunction(15)
+print *, "function should be:", sourcePoints(10,1)**2*sourcePoints(10,2), sourcePoints(10,3), sourcePoints(10,1)
+print *, "function value dot normal:", sum(functionValue*unitVectors(10,:))
+print *, "function dot normal should be:", testFunction(10)
- end subroutine test_ibInterp_computeVectorNoSlipRBFCoefficients
-!
-! subroutine test_ibInterp_computeVectorNoSlipPlanarRBFCoefficients(pointCount, &
-! sourcePoints, unitVectors, destinationPoint, &
-! alpha, planeBasisVectors, coefficients, dminfo)
-!
-! 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
-! type (dm_info), intent(in) :: dminfo
-!
-! 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
-! real(kind=RKIND), dimension(:,:), pointer :: rhs, coeffs
-! integer, dimension(:), pointer :: pivotIndices
-!
-! matrixSize = pointCount+2 ! space for constant vector in plane
-!
-! allocate(matrix(matrixSize,matrixSize))
-! allocate(rhs(matrixSize,2))
-! allocate(coeffs(matrixSize,2))
-! allocate(pivotIndices(matrixSize))
-!
-! matrix = 0.0
-! rhs = 0.0
-! coeffs = 0.0
-!
-! do i = 1, pointCount
-! planarSourcePoints(i,1) = sum(sourcePoints(i,:)*planeBasisVectors(1,:))
-! planarSourcePoints(i,2) = sum(sourcePoints(i,:)*planeBasisVectors(2,:))
-! planarUnitVectors(i,1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
-! planarUnitVectors(i,2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
-! end do
-! planarDestinationPoint(1) = sum(destinationPoint*planeBasisVectors(1,:))
-! planarDestinationPoint(2) = sum(destinationPoint*planeBasisVectors(2,:))
-! call setUpVectorNoSlipRBFMatrixAndRHS(pointCount, 2, &
-! planarSourcePoints, planarUnitVectors, planarDestinationPoint, &
-! alpha, matrix, rhs)
-!
-! do i = 1, pointCount
-! matrix(i,pointCount+1) = sum(unitVectors(i,:)*planeBasisVectors(1,:))
-! matrix(i,pointCount+2) = sum(unitVectors(i,:)*planeBasisVectors(2,:))
-! matrix(pointCount+1:pointCount+2,i) = matrix(i,pointCount+1:pointCount+2)
-! end do
-! rhs(pointCount+1,:) = planeBasisVectors(1,:)
-! rhs(pointCount+2,:) = planeBasisVectors(2,:)
-!
-! ! solve each linear system
-! call LEGS(matrix, matrixSize, rhs(:,1), coeffs(:,1), pivotIndices)
-! call LEGS(matrix, matrixSize, rhs(:,2), coeffs(:,2), pivotIndices)
-!
-! coefficients(:,1) = planeBasisVectors(1,1)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,1)*coeffs(1:pointCount,2)
-! coefficients(:,2) = planeBasisVectors(1,2)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,2)*coeffs(1:pointCount,2)
-! coefficients(:,3) = planeBasisVectors(1,3)*coeffs(1:pointCount,1) &
-! + planeBasisVectors(2,3)*coeffs(1:pointCount,2)
-!
-! deallocate(matrix)
-! deallocate(rhs)
-! deallocate(coeffs)
-! deallocate(pivotIndices)
-!
-! end subroutine test_ibInterp_computeVectorNoSlipPlanarRBFCoefficients
-!
+ end subroutine test_ibInterp_computeVectorDirichletRBFCoefficients
+
+ subroutine test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
+ integer, parameter :: pointCount = 14
+ real(kind=RKIND), dimension(pointCount,3) :: sourcePoints
+ real(kind=RKIND), dimension(pointCount,3) :: unitVectors
+ real(kind=RKIND), dimension(3) :: destinationPoint, functionValue
+ real(kind=RKIND) :: alpha
+ real(kind=RKIND), dimension(2,3) :: planeBasisVectors
+ real(kind=RKIND), dimension(pointCount, 3) :: coefficients
+
+ real(kind=RKIND), dimension(pointCount) :: testFunction
+
+ alpha = 1.35
+
+ planeBasisVectors(1,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ planeBasisVectors(2,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+
+
+ sourcePoints(1,:) = (/-0.408248290463863, -0.502982635954617, 0.911230926418480/)
+ sourcePoints(2,:) = (/-0.816496580927726, 0.0546948998705890, 0.761801681057137/)
+ sourcePoints(3,:) = (/0.408248290463863, -0.911230926418480, 0.502982635954617/)
+ sourcePoints(4,:) = (/0.0, -0.353553390593274, 0.353553390593274/)
+ sourcePoints(5,:) = (/0.816496580927726, -0.761801681057137, -0.0546948998705890/)
+ sourcePoints(6,:) = (/-0.408248290463863, 0.204124145231932, 0.204124145231932/)
+ sourcePoints(7,:) = (/-0.816496580927726, 0.761801681057137, 0.0546948998705890/)
+ sourcePoints(8,:) = (/0.408248290463863, -0.204124145231932, -0.204124145231932/)
+ sourcePoints(9,:) = (/0.0, 0.353553390593274, -0.353553390593274/)
+ sourcePoints(10,:) = (/0.734846922834954, -0.155291427061512, -0.579555495773441/)
+ sourcePoints(11,:) = (/-0.408248290463863, 0.911230926418480, -0.502982635954617/)
+ sourcePoints(12,:) = (/0.244948974278318, 0.301789581572770, -0.546738555851088/)
+ sourcePoints(13,:) = (/0.734846922834954, -0.155291427061512, -0.579555495773441/)
+ sourcePoints(14,:) = (/0.244948974278318, 0.301789581572770, -0.546738555851088/)
+
+
+ unitVectors(1,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ unitVectors(2,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ unitVectors(3,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ unitVectors(4,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ unitVectors(5,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ unitVectors(6,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ unitVectors(7,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ unitVectors(8,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ unitVectors(9,:) = (/0.0, 0.707106781186548, -0.707106781186548/)
+ unitVectors(10,:) = (/0.365148371670111, 0.449881346198621, -0.815029717868732/)
+ unitVectors(11,:) = (/0.816496580927726, -0.408248290463863, -0.408248290463863/)
+ unitVectors(12,:) = (/0.365148371670111, 0.449881346198621, -0.815029717868732/)
+ unitVectors(13,:) = (/0.730296743340221, -0.681376137686949, -0.0489206056532726/)
+ unitVectors(14,:) = (/0.730296743340221, -0.681376137686949, -0.0489206056532726/)
+
+ testFunction = (/-0.250000000000000, -0.500000000000000, -0.250000000000000, -0.500000000000000, -0.500000000000000, 0.0, 0.500000000000000, 0.0, 0.500000000000000, 0.377001061006465, 0.250000000000000, 0.560805848756947, 0.0831817287629922, -0.220029088985979/)
+
+ destinationPoint = sourcePoints(7,:)
+
+ call ibInterp_computeVectorDirichletPlanarRBFCoefficients(pointCount, &
+ sourcePoints, unitVectors, destinationPoint, &
+ alpha, planeBasisVectors, coefficients)
+
+print *, "coefficients", coefficients
+
+ functionValue(1) = sum(coefficients(:,1)*testFunction)
+ functionValue(2) = sum(coefficients(:,2)*testFunction)
+ functionValue(3) = sum(coefficients(:,3)*testFunction)
+
+print *, "function value:", functionValue
+print *, "function should be: 0.408248290463863 0.149429245361343 -0.557677535825206"
+print *, "function value dot normal:", sum(functionValue*unitVectors(7,:))
+print *, "function dot normal should be:", testFunction(7)
+
+ end subroutine test_ibInterp_computeVectorDirichletPlanarRBFCoefficients
+
! subroutine test_ibInterp_computeVectorFreeSlipRBFCoefficients(pointCount, &
! sourcePoints, isTangentToInterface, normalVectorIndex, unitVectors, destinationPoint, &
! alpha, coefficients, dminfo)
Modified: branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F
===================================================================
--- branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-04-28 23:25:27 UTC (rev 222)
+++ branches/ocean_projects/IBinterp/mpas/src/operators/module_vector_reconstruction.F 2010-04-29 04:02:42 UTC (rev 223)
@@ -3,10 +3,11 @@
use grid_types
use configure
use constants
+ use IB_interpolation
implicit none
- public :: init_reconstruct, reconstruct
+ public :: init_reconstruct, reconstruct, new_init_reconstruct
contains
@@ -102,6 +103,7 @@
normals(:,j) = X1(:) - X2(:)
endif
+
call unit_vector_in_R3(normals(:,j))
xLoc(:,j,iCell) = 0.5*(X2(:) + X1(:))
@@ -128,6 +130,7 @@
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
@@ -146,8 +149,7 @@
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))
@@ -176,6 +178,7 @@
+ rbfValue * normals(:,i) * invMatrix(i,j)
enddo
enddo
+
! polynomial parts
do j=1,npoints
coeffs_reconstruct(:,j,iCell) = coeffs_reconstruct(:,j,iCell) &
@@ -195,65 +198,161 @@
end subroutine init_reconstruct
+ subroutine new_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
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine reconstruct(s, grid)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! 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
- implicit none
+ type (grid_meta), intent(inout) :: grid
- type (grid_state), intent(inout) :: s
- type (grid_meta), intent(in) :: grid
+ ! temporary arrays needed in the (to be constructed) init procedure
+ integer :: nCellsSolve
+ integer, dimension(:,:), pointer :: edgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
+ integer :: i, iCell, pointCount, maxEdgeCount
+ real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
+ real (kind=RKIND) :: r, cellCenter(3), alpha, tangentPlane(2,3)
+ real (kind=RKIND), allocatable, dimension(:,:) :: edgeOnCellLocations, edgeOnCellNormals, &
+ coeffs
- ! 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 :: u
- real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+ real(kind=RKIND), dimension(:,:), pointer :: edgeNormalVectors, edgeLocations
+ real(kind=RKIND), dimension(:,:,:), pointer :: cellTangentPlane
- real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
- ! stored arrays used during compute procedure
- coeffs_reconstruct => grid % coeffs_reconstruct % array
+ !========================================================
+ ! arrays filled and saved during init procedure
+ !========================================================
+ coeffs_reconstruct => grid % coeffs_reconstruct % array
- ! temporary variables
- edgesOnCell => grid % edgesOnCell % array
- nEdgesOnCell=> grid % nEdgesOnCell % array
- nCellsSolve = grid % nCellsSolve
- u => s % u % array
- uReconstructX => s % uReconstructX % array
- uReconstructY => s % uReconstructY % array
- uReconstructZ => s % uReconstructZ % array
+ !========================================================
+ ! temporary variables needed for init procedure
+ !========================================================
+ xCell => grid % xCell % array
+ yCell => grid % yCell % array
+ zCell => grid % zCell % array
+ edgesOnCell => grid % edgesOnCell % array
+ nEdgesOnCell=> grid % nEdgesOnCell % array
+ nCellsSolve = grid % nCellsSolve
+ edgeNormalVectors => grid % edgeNormalVectors % array
+ edgeLocations => grid % edgeLocations % array
+ cellTangentPlane => grid % cellTangentPlane % array
- ! 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
+ ! init arrays
+ coeffs_reconstruct = 0.0
- end subroutine reconstruct
+ 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
+ edgeOnCellLocations(i,:) = edgeLocations(:, edgesOnCell(i,iCell))
+ edgeOnCellNormals(i,:) = edgeNormalVectors(:, edgesOnCell(i,iCell))
+ 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 ibInterp_computeVectorDirichletPlanarRBFCoefficients(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 new_init_reconstruct
+
+ subroutine reconstruct(state, grid)
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! 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 (grid_state), intent(inout) :: state
+ type (grid_meta), intent(in) :: grid
+
+ ! 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 :: u
+ real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: coeffs_reconstruct
+
+ ! 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
+ u => state % u % array
+ uReconstructX => state % uReconstructX % array
+ uReconstructY => state % uReconstructY % array
+ uReconstructZ => state % uReconstructZ % array
+
+ ! 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
+
+ end subroutine reconstruct
+
+
subroutine get_distance(x1,x2,xout)
implicit none
real(kind=RKIND), intent(in) :: x1(3), x2(3)
</font>
</pre>