most recent fortran development sync.

development/BIVARParticlesInterface.f90

@@ -9,53 +9,161 @@ module BIVARInterfaceModule
 implicit none
 
 private
-public BIVARInterface, New, Delete, Copy
-public InterpolateScalarField, InterpolateVectorField
-public EstimateDerivatives
+public BIVARInterface, New, Delete
+public InterpolateScalar, InterpolateVector
 
 !
 !----------------
 ! Types and module constants
 !----------------
 !
 type BIVARInterface
- integer(kint) :: nSourcePoints = 0
- integer(kint) :: nTri = 0
- integer(kint), pointer :: triVerts(:) => null()
- integer(kint) :: nBoundaryTriEdges = 0
- integer(kint), pointer :: boundaryEdgesAndTri(:) => null()
- real(kreal), pointer :: scalarPartials(:) => null()
+ real(kreal), dimension(:), pointer :: estPartials1 => null()
+ real(kreal), dimension(:), pointer :: estPartials2 => null()
 end type
 
-interface New
- module procedure newPrivate
-end interface
+type PlaneTri
+ integer(kint) :: nTri
+ integer(kint), pointer, dimension(:) :: vertIndices => null()
+ integer(kint) :: nBorderSegments
+ integer(kint), pointer, dimension(:) :: borderIndices => null()
+ integer(kint), pointer, dimension(:) :: integerWork => null()
+ real(kreal), pointer, dimension(:) :: realWork => null()
+ contains
+ final :: deleteTri
+end type
 
-interface Copy
- module procedure copyPrivate
+interface New
+ module procedure newPrivateField
+ module procedure newPrivateLagParam
+ module procedure newTri
 end interface
 
 interface Delete
  module procedure deletePrivate
+ module procedure deleteTri
 end interface
 
 contains
 
-subroutine newPrivate( self, sourceParticles, sourceField )
+subroutine newTri( delTri, aParticles )
+ type(PlaneTri), intent(out) :: delTri
+ type(Particles), intent(in) :: aParticles
+ !
+ integer(kint), allocatable, dimension(:) :: iwp
+
+ allocate(delTri%vertIndices( 6 * aParticles%N - 15) )
+ allocate(delTri%borderIndices(6* aParticles%N) )
+ allocate(delTri%integerWork( 18 * aParticles%N))
+ allocate(delTri%realWork(aParticles%N))
+
+ allocate(iwp( aParticles%N))
+
+ call IDTANG( aParticles%N, aParticles%x, aParticles%y, delTri%nTri, delTri%vertIndices, delTri%nBorderSegments, &
+ delTri%borderIndices, delTri%integerWork, iwp, delTri%realWork)
+
+ deallocate(iwp)
+end subroutine
+
+subroutine deleteTri(delTri)
+ type(PlaneTri), intent(inout) :: delTri
+ if ( associated(delTri%vertIndices)) deallocate(delTri%vertIndices)
+ if ( associated(delTri%borderIndices)) deallocate(delTri%borderIndices)
+ if ( associated(delTri%integerWork)) deallocate(delTri%integerWork)
+ if ( associated(delTri%realWork)) deallocate(delTri%realWork)
+end subroutine
+
+subroutine newPrivateField( self, delTri, sourceParticles, sourceField )
  type(BIVARInterface), intent(out) :: self
+ type(PlaneTri), intent(in) :: delTri
  type(Particles), intent(in) :: sourceParticles
  type(Field), intent(in) :: sourceField
+ !
+ real(kreal), allocatable, dimension(:) :: realWork
+
+ allocate(self%estPartials1( 5 * sourceParticles%N ))
+ if ( sourceField%nDim == 2 ) allocate(self%estPartials2( 5 * sourceParticles%N))
+
+ allocate(realWork(sourceParticles%N))
 
+ if ( sourceField%nDim == 1 ) then
+ call IDPDRV( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%scalar, &
+ delTri%nTri, delTri%vertIndices, self%estPartials1, realWork)
+ elseif ( sourceField%nDim == 2 ) then
+ call IDPDRV( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%xComp, &
+ delTri%nTri, delTri%vertIndices, self%estPartials1, realWork)
+ call IDPDRV( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%yComp, &
+ delTri%nTri, delTri%vertIndices, self%estPartials1, realWork) 
+ endif
+
+ deallocate(realWork)
 end subroutine
 
-subroutine InterpolateScalarField( aParticles, scalarField, xOut, yOut, interpOut )
- type(Particles), intent(in) :: aParticles
- type(Field), intent(in) :: scalarField
- real(kreal), intent(in) :: xOut(:)
- real(kreal), intent(in) :: yOut(:)
- real(kreal), intent(inout) :: interpOut(:)
+subroutine newPrivateLagParam( self, delTri, sourceParticles )
+ type(BIVARInterface), intent(out) :: self
+ type(PlaneTri), intent(in) :: delTri
+ type(Particles), intent(in) :: sourceParticles
+ !
+ real(kreal), allocatable, dimension(:) :: realWork
+
+ allocate(self%estPartials1( 5 * sourceParticles%N ))
+ allocate(self%estPartials2( 5 * sourceParticles%N ))
 
+ allocate(realWork(sourceParticles%N))
 
+ call IDPDRV( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceParticles%x0, &
+ delTri%nTri, delTri%vertIndices, self%estPartials1, realWork)
+ call IDPDRV( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceParticles%y0, &
+ delTri%nTri, delTri%vertIndices, self%estPartials1, realWork) 
+
+ deallocate(realWork)
 end subroutine
 
+subroutine deletePrivate(self)
+ type(BIVARInterface), intent(inout) :: self
+ if ( associated(self%estPartials1)) deallocate(self%estPartials1)
+ if ( associated(self%estPartials2)) deallocate(self%estPartials2)
+end subroutine
+
+function InterpolateScalar( self, delTri, sourceParticles, sourceField, xOut, yOut )
+ real(kreal) :: InterpolateScalar 
+ type(BIVARInterface), intent(in) :: self
+ type(PlaneTri), intent(in) :: delTri
+ type(Particles), intent(in) :: sourceParticles
+ type(Field), intent(in) :: sourceField
+ real(kreal), intent(in) :: xOut
+ real(kreal), intent(in) :: yOut
+ !
+ integer(kint) :: inTri
+
+ call IDLCTN( sourceParticles%N, sourceParticles%x, sourceParticles%y, delTri%nTri, delTri%vertIndices, &
+ delTri%nBorderSegments, delTri%borderIndices, xOut, yOut, inTri, delTri%integerWork, delTri%realWork)
+
+ call IDPTIP( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%scalar, &
+ delTri%nTri, delTri%vertIndices, delTri%nBorderSegments, delTri%borderIndices, self%estPartials1, inTri, &
+ xOut, yOut, InterpolateScalar)
+end function
+
+function InterpolateVector( self, delTri, sourceParticles, sourceField, xOut, yOut )
+ real(kreal), dimension(2) :: InterpolateVector 
+ type(BIVARInterface), intent(in) :: self
+ type(PlaneTri), intent(in) :: delTri
+ type(Particles), intent(in) :: sourceParticles
+ type(Field), intent(in) :: sourceField
+ real(kreal), intent(in) :: xOut
+ real(kreal), intent(in) :: yOut
+ !
+ integer(kint) :: inTri
+
+ call IDLCTN( sourceParticles%N, sourceParticles%x, sourceParticles%y, delTri%nTri, delTri%vertIndices, &
+ delTri%nBorderSegments, delTri%borderIndices, xOut, yOut, inTri, delTri%integerWork, delTri%realWork)
+
+ call IDPTIP( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%xComp, &
+ delTri%nTri, delTri%vertIndices, delTri%nBorderSegments, delTri%borderIndices, self%estPartials1, inTri, &
+ xOut, yOut, InterpolateVector(1))
+ call IDPTIP( sourceParticles%N, sourceParticles%x, sourceParticles%y, sourceField%yComp, &
+ delTri%nTri, delTri%vertIndices, delTri%nBorderSegments, delTri%borderIndices, self%estPartials2, inTri, &
+ xOut, yOut, InterpolateVector(2))
+end function 
+
 end module

development/CMakeLists.txt

@@ -41,7 +41,8 @@ add_library( lpmFortran NumberKinds.f90 OutputWriter.f90 Logger.f90 IntegerList.
  SphereGeometry.f90 PlaneGeometry.f90 
  Particles.f90 Field.f90
  Edges.f90 Faces.f90 PolyMesh2d.f90 
- ScatteredTriInterp.f90 bivar.f90 PSE.f90 PSEDirectSum.f90 MLSQ.f90 TrivariateQuadraticMLSQ.f90
+ ScatteredTriInterp.f90 bivar.f90 BIVARParticlesInterface.f90
+ PSE.f90 PSEDirectSum.f90 MLSQ.f90 TrivariateQuadraticMLSQ.f90
  stripack.f ssrfpack.f ssrfpackInterface.f90
  MPISetup.f90 PlanarSWE.f90 SWEPlaneRK4.f90 
  Refinement.f90 SphereSWE.f90 SphereSWESolver.f90
@@ -111,6 +112,9 @@ TARGET_LINK_LIBRARIES( ssrfpackTest.exe lpmFortran )
 ADD_EXECUTABLE( sphereTrivarMLSQ.exe TrivarMLSQSphereTest.f90)
 TARGET_LINK_LIBRARIES( sphereTrivarMLSQ.exe lpmFortran ${LAPACK_LIBRARIES})
 
+ADD_EXECUTABLE( spherePSEInterpTest.exe spherePSEInterpTest.f90 )
+TARGET_LINK_LIBRARIES( spherePSEInterpTest.exe lpmFortran)
+
 # Unit Tests
 ENABLE_TESTING()
 ADD_TEST( particlesAndFieldTest ${PROJECT_BINARY_DIR}/particlesAndFieldTest.exe)

development/SSRFPACKTest.f90

@@ -104,12 +104,12 @@ program ssrfpacktest
 
 if ( meshSeed == ICOS_TRI_SPHERE_SEED ) then
  write(logString,'(A,I1,A)') "building icos tri mesh to initNest = ", initNest, "..."
- write(paraviewFilename,'(A,I1,A)') "icosTriPSE_", initNest, ".vtk"
- write(matlabFilename,'(A,I1,A)') "icosTriPSE_", initNest, ".m"
+ write(paraviewFilename,'(A,I1,A)') "icosTriSSRFPACK_", initNest, ".vtk"
+ write(matlabFilename,'(A,I1,A)') "icosTriSSRFPACK_", initNest, ".m"
 else
  write(logString,'(A,I1,A)') "building cubed sphere mesh to initNest = ", initNest, "..."
- write(paraviewFilename,'(A,I1,A)') "cubedSpherePSE_", initNest, ".vtk"
- write(matlabFilename,'(A,I1,A)') "cubedSpherePSE_", initNest, ".m"
+ write(paraviewFilename,'(A,I1,A)') "cubedSphereSSRFPACK_", initNest, ".vtk"
+ write(matlabFilename,'(A,I1,A)') "cubedSphereSSRFPACK_", initNest, ".m"
 endif
 call LogMessage(exeLog, TRACE_LOGGING_LEVEL, logKey, logString )
 
@@ -174,7 +174,7 @@ program ssrfpacktest
 do j = unifGridMPI%indexStart(procRank), unifGridMPI%indexEnd(procRank)
  do i = 1, 181
  interp(i,j) = InterpolateScalar( lons(j), lats(i), ssrfConst, sphere, delTri, constantScalar)
- harmInterp(i,j) = InterpolateScalar(lons(j), lats(i), ssrfHarm, sphere, delTri, constantScalar)
+ harmInterp(i,j) = InterpolateScalar(lons(j), lats(i), ssrfHarm, sphere, delTri, sphHarm)
  enddo
 enddo