Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_test) Browse

fennecs/src/splinebound_mod.f90
77b4b66b5239geom_test

Last Change
Raw File

splinebound_mod.f90
View Options

MODULE splinebound
USE constants
USE bsplines
USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints
Use forSISLdata
IMPLICIT NONE
INTEGER(SELECTED_INT_KIND(4)), PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2
TYPE spline_boundary
! all curves assume right handedness to set which side of the curve is inside or outside
type(SISLCurve):: curve
type(SISLCurve):: bdcurve
Real(kind=db):: Dirichlet_val !< Value for the dirichlet boundary condition created by this boundary
Real(kind=db):: epsge=1.0e-5 !< geometric resolution used for calculating distances
Real(kind=db):: epsce=1.0e-9 !< value of weight below which it is 0
INTEGER(kind(bd)):: type=bd_Dirichletconst !< type of boundary conditions
END TYPE spline_boundary
type spline_domain
integer:: nbsplines !< number of spline boundaries in the domain
type(spline_boundary), allocatable:: boundaries(:) !< List of boundaries in the domain
Real(kind=db):: dist_extent=0.1 !< distance used for the merging with the plateau function for the weight
type(cellkind), ALLOCATABLE:: cellk(:,:) !< Precomputed parameters at each cell for faster weight computation
type(spline2d), pointer:: splrz => null() !< Pointer to the main spline grid used for the FEM solver
Integer:: nb1 !< Number of grid points in the 1st dimension
Integer:: nb2 !< Number of grid points in the 2nd dimension
real(kind=db), ALLOCATABLE:: x1(:) !< Grid points in first direction for weight interpolation
real(kind=db), ALLOCATABLE:: x2(:) !< Grid points in 2nd direction for weight interpolation
real(kind=db), ALLOCATABLE:: dx1(:) !< inverse cell width in first direction for weight interpolation
real(kind=db), ALLOCATABLE:: dx2(:) !< inverse cell width in 2nd direction for weight interpolation
type(SISLsurf):: Dirdomweight
type(SISLsurf):: totdomweight
end type spline_domain
type cellkind
integer:: spldirkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate
integer:: spltotkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate
integer:: linkedboundaries(2)=-1 !< stores the spline curve indices in the spline_domain of the spline boundaries that are at a distance lower than dist_extent
integer:: uniontype=0 !< defines if the two boundaries are linked by 1: union, -1: intersection or 0, no joint
integer:: leftknot(2)=0 !< for each linkedboundary saves a lower initial guess for calculating the distance
real(kind=db):: lguess(2)=-1 !< Spline parameter left limit as start guess
real(kind=db):: rguess(2)=-1 !< Spline parameter right limit as start guess
end type cellkind
CONTAINS
subroutine read_splinebound(Fileid, spldom, splrz, rnorm, Phinorm, Potinn,Potout)
use mpi
Integer:: Fileid
type(spline_domain):: spldom
type(spline2d):: splrz
real(kind=db):: rnorm, phinorm
Integer:: nbsplines, istat, mpirank, ierr
real(kind=db):: dist_extent, Potinn, Potout, epsge, epsce
Character(len=128):: h5fname="", line
real(kind=db) :: Dvals(30)=0
integer:: nbpoints
real(kind=db),allocatable:: cpoints(:,:)
real(kind=db):: dz, dr, ctr, Lz, ra, rb
integer:: degree=4, i
namelist /spldomain/ nbsplines, dist_extent, h5fname, dr, Lz, ra, rb, epsge, epsce, Dvals
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
REWIND(fileid)
READ(fileid,spldomain, iostat=istat)
if (istat.gt.0) then
if(mpirank .eq. 0) then
backspace(fileid)
read(fileid,fmt='(A)') line
write(*,'(A)') &
'Invalid line in geomparams: '//trim(line)
end if
call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
stop
end if
if(mpirank .eq. 0) WRITE(*, spldomain)
Dvals=Dvals/phinorm
dist_extent=dist_extent/rnorm
if (.not. trim(h5fname)=='' ) then
call setspline_domain(spldom, splrz, dist_extent, 0)
call splinebound_readh5domain(h5fname,spldom, rnorm, phinorm)
call classifycells(spldom)
do i=1,spldom%nbsplines
spldom%boundaries(i)%Dirichlet_val=Dvals(i)
end do
return
end if
nbpoints=600
allocate(cpoints(2,nbpoints))
nbsplines=2
call setspline_domain(spldom, splrz, dist_extent, nbsplines)
dz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/(nbpoints-1)
Lz=Lz/rnorm
ctr=(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0))/2
do i=1,nbpoints
cpoints(1,i)=splrz%sp1%knots(0)+(i-1)*dz
end do
cpoints(2,:)=ra/rnorm
call setspline_boundary(spldom%boundaries(1), cpoints, degree, Potinn, epsge, epsce)
!spldom%boundaries(1)%type=bd_Neumann
do i=1,nbpoints
cpoints(1,i)=splrz%sp1%knots(0)+(nbpoints-i)*dz
end do
dr=dr/rnorm
cpoints(2,:)=rb/rnorm
do i=1,nbpoints
if(((cpoints(1,i)-ctr)/Lz)**2.le.1) then
cpoints(2,i)=rb/rnorm+dr*sqrt(1-(cpoints(1,i)-ctr)**2/Lz**2)
endif
end do
call setspline_boundary(spldom%boundaries(2), cpoints, degree, Potout, epsge, epsce)
do i=1,spldom%nbsplines
spldom%boundaries(i)%Dirichlet_val=Dvals(i)
end do
call classifycells(spldom)
end subroutine
Subroutine splinebound_diag(File_handle, str, spldom)
Use mpi
Use futils
Use basic, ONLY: rnorm, phinorm
Integer:: File_handle
type(spline_domain):: spldom
Character(len=*):: str
CHARACTER(len=128):: grpname
Integer:: ierr, mpirank, i
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
IF(mpirank .eq. 0) THEN
Write(grpname,'(a,a)') trim(str),"/geometry_spl"
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "dist_extent",spldom%dist_extent)
Call attach(File_handle, trim(grpname), "nbsplines", spldom%nbsplines)
do i=1,spldom%nbsplines
Write(grpname,'(a,a,i2.2)') trim(str),"/geometry_spl/",i
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "Dirichlet_val", spldom%boundaries(i)%Dirichlet_val*phinorm)
Call attach(File_handle, trim(grpname), "order", spldom%boundaries(i)%curve%ik)
Call attach(File_handle, trim(grpname), "kind", spldom%boundaries(i)%curve%ikind)
Call attach(File_handle, trim(grpname), "dim", spldom%boundaries(i)%curve%idim)
CALL putarr(File_handle, TRIM(grpname)//"/pos", spldom%boundaries(i)%curve%ecoef*rnorm)
CALL putarr(File_handle, TRIM(grpname)//"/knots", spldom%boundaries(i)%curve%et)
end do
END IF
End subroutine splinebound_diag
subroutine splinebound_readh5domain(filename, spldom, rnorm, phinorm)
use futils
use forSISL
implicit none
Character(len=*),intent(in) :: filename
type(spline_domain),intent(inout) :: spldom
integer:: h5id, i
real(kind=db):: rnorm, phinorm
CHARACTER(len=128):: grpname
logical:: periodic
integer:: order, dim, bdtype
INTEGER:: posrank, posdim(2)
real(kind=db):: Dval, epsge, epsce
real(kind=db),allocatable:: points(:,:)
call openf(filename, h5id,'r','d')
call getatt(h5id, '/geometry_spl/','nbsplines', spldom%nbsplines)
!call getatt(h5id, '/geometry_spl/','dist_extent', spldom%dist_extent)
!spldom%dist_extent=spldom%dist_extent/rnorm
! prepare memory
if (allocated(spldom%boundaries)) then
do i=1,size(spldom%boundaries,1)
call free_bsplinecurve(spldom%boundaries(i))
end do
DEALLOCATE(spldom%boundaries)
end if
allocate(spldom%boundaries(spldom%nbsplines))
! Read each boundary curve individually
do i=1,spldom%nbsplines
Write(grpname,'(a,i2.2)') "/geometry_spl/",i
If(.not. isgroup(h5id, trim(grpname))) THEN
Write(*,*) "Error the geometry definition file is invalid"
END IF
Call getatt(h5id, trim(grpname), "Dirichlet_val", Dval)
Call getatt(h5id, trim(grpname), "epsge", epsge)
Call getatt(h5id, trim(grpname), "epsce", epsce)
Call getatt(h5id, trim(grpname), "order", order)
Call getatt(h5id, trim(grpname), "dim", dim)
Call getatt(h5id, trim(grpname), "periodic", periodic)
CALL getdims(h5id, TRIM(grpname)//"/pos", posrank, posdim)
allocate(points(posdim(1),posdim(2)))
CALL getarr(h5id, TRIM(grpname)//"/pos", points)
points=points/rnorm
Call setspline_boundary(spldom%boundaries(i),transpose(points), order-1, Dval/phinorm, epsge,epsce, periodic)
bdtype=bd
Call getatt(h5id, trim(grpname), "type", bdtype)
if(bdtype .ne. -1) spldom%boundaries(i)%type=bdtype
deallocate(points)
end do
call closef(h5id)
end subroutine splinebound_readh5domain
subroutine setspline_domain(spldom,splrz,dist_extent, nb_splines)
type(spline_domain):: spldom
type(spline2d), TARGET:: splrz
real(kind=db):: dist_extent
integer:: nb_splines, nb1, nb2
! Store the grid parameters to speed-up calculations
nb1=splrz%sp1%nints
nb2=splrz%sp2%nints
spldom%nb1=nb1
spldom%nb2=nb2
spldom%splrz=>splrz
allocate(spldom%cellk(0:nb1-1,0:nb2-1))
allocate(spldom%x1(0:nb1))
allocate(spldom%x2(0:nb2))
allocate(spldom%dx1(0:nb1-1))
allocate(spldom%dx2(0:nb2-1))
spldom%x1(0:)=splrz%sp1%knots(0:nb1)
spldom%x2(0:)=splrz%sp2%knots(0:nb2)
spldom%dx1(0:)=1/(spldom%x1(1:nb1)-spldom%x1(0:nb1-1))
spldom%dx2(0:)=1/(spldom%x2(1:nb2)-spldom%x2(0:nb2-1))
!Prepare structures to host singular spline boundaries
spldom%nbsplines=nb_splines
if(spldom%nbsplines.gt. 0) allocate(spldom%boundaries(nb_splines))
spldom%dist_extent=dist_extent
end subroutine setspline_domain
subroutine setspline_boundary(b_curve, cpoints, degree, D_val, epsge, epsce, periodic)
Use bsplines
use mpi
type(spline_boundary):: b_curve
Real(kind=db):: cpoints(:,:)
Real(REAL64),ALLOCATABLE:: knots(:),points(:)
integer:: degree
logical, optional:: periodic
Integer:: order, ierr
Real(kind=db):: D_val
Real(kind=db),OPTIONAL :: epsge, epsce
Integer:: nbpoints,i, dim, kind, copy
logical:: period
period=.false.
if(present(periodic))period=periodic
nbpoints= size(cpoints,2)
dim=size(cpoints,1)
order=degree+1
if(nbpoints .lt. order) then
WRITE(*,'(a,i3,a,i5)') "Error: the number of points", nbpoints, " is insuficient for the required order ", order
CALL mpi_finalize(ierr)
call EXIT(-1)
end if
kind=1 ! polynomial b-spline curve
copy = 1 ! copy input array
Allocate(knots(1:nbpoints+order))
DO i=order,nbpoints+1
knots(i)=real(i-order)/real(nbpoints+1-order)
END DO
knots(1:order)=knots(order)
knots(nbpoints+1:nbpoints+order)=knots(nbpoints+1)
allocate(points(dim*nbpoints))
points=reshape(cpoints,(/dim*nbpoints/))
CALL newCurve(nbpoints, order, knots, points, kind, dim, copy, b_curve%curve)
b_curve%Dirichlet_val=D_val
if(present(epsge)) b_curve%epsge=epsge
if(present(epsce)) b_curve%epsce=epsce
end subroutine setspline_boundary
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Calculates the geometric weight from a given spline domain
!> @param[in] spldom spline domain containing the information on the boundary conditions
!> @param[in] x1(:) array of axial positions where the weights are evaluated
!> @param[in] x2(:) array of radial positions where the weights are evaluated
!> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative
!---------------------------------------------------------------------------
SUBROUTINE spline_w(spldom,x1,x2,w)
use forSISL,ONLY: s1424
type(spline_domain):: spldom
Real(kind=db), INTENT(IN):: x2(:),x1(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Integer:: k,sstat,l
type(cellkind):: cellk
REAL(real64):: wtmp(4), point(2)
Integer,allocatable::i(:),j(:)
allocate(i(size(x2,1)),j(size(x2,1)))
call getindex(x1, x2, spldom, i, j)
Do k=1,size(x2,1)
cellk=spldom%cellk(i(k),j(k))
SELECT CASE (cellk%spldirkind)
CASE(-1)
w(k,:)=0
w(k,0)=-1
CYCLE
CASE(1)
w(k,:)=0
w(k,0)=1
CYCLE
CASE DEFAULT
!call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1))
point=(/x1(k),x2(k)/)
call s1424(spldom%Dirdomweight,1,1,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat)
if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight"
if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight"
w(k,:)=wtmp(1:size(w,2))
END SELECT
end DO
End SUBROUTINE spline_w
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Calculates the geometric weight from a given spline domain
!> @param[in] spldom spline domain containing the information on the boundary conditions
!> @param[in] x1(:) array of axial positions where the weights are evaluated
!> @param[in] x2(:) array of radial positions where the weights are evaluated
!> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative
!---------------------------------------------------------------------------
SUBROUTINE spline_wtot(spldom,x1,x2,w)
use forSISL,ONLY: s1424
type(spline_domain):: spldom
Real(kind=db), INTENT(IN):: x2(:),x1(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Integer:: k,sstat,l
type(cellkind):: cellk
REAL(real64):: wtmp(4), point(2)
Integer,allocatable::i(:),j(:)
allocate(i(size(x2,1)),j(size(x2,1)))
call getindex(x1, x2, spldom, i, j)
Do k=1,size(x2,1)
cellk=spldom%cellk(i(k),j(k))
SELECT CASE (cellk%spltotkind)
CASE(-1)
w(k,:)=0
w(k,0)=-1
CYCLE
CASE(1)
w(k,:)=0
w(k,0)=1
CYCLE
CASE DEFAULT
!call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1))
point=(/x1(k),x2(k)/)
if(size(w,2).gt.1) then
call s1424(spldom%totdomweight,1,1,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat)
else
call s1424(spldom%totdomweight,0,0,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat)
end if
if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splinedomweight"
if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splinedomweight"
w(k,:)=wtmp(1:size(w,2))
END SELECT
end DO
End SUBROUTINE spline_wtot
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Calculates the interpolation in the domain of the Dirichlet boundary conditions from a given spline domain
!> @param[in] spldom spline domain containing the information on the boundary conditions
!> @param[in] z(:) array of axial positions where the weights are evaluated
!> @param[in] r(:) array of radial positions where the weights are evaluated
!> @param[out] g(:,0:) matrix of boundary interpolations g with first index corresponding to the position and second index to the derivative
!---------------------------------------------------------------------------
SUBROUTINE spline_g(spldom,x1,x2,g,w)
use forSISL,ONLY: s1424
type(spline_domain):: spldom
Real(kind=db), INTENT(IN):: x2(:),x1(:)
Real(kind=db), INTENT(OUT):: g(:,0:)
Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:)
REAL(real64),allocatable:: gtmp(:,:)
Integer:: k,sstat,sizew
Integer,allocatable::i(:),j(:)
type(cellkind):: cellk
allocate(gtmp(size(x2,1),0:size(g,2)-1))
allocate(i(size(x2,1)),j(size(x2,1)))
call getindex(x1, x2, spldom, i, j)
if(present(w)) then
gtmp=w
else
Do k=1,size(x2,1)
cellk=spldom%cellk(i(k),j(k))
call s1424(spldom%Dirdomweight,1,1,(/x1(k),x2(k)/),cellk%leftknot(1),cellk%leftknot(2),gtmp(k,:),sstat)
if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight"
if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight"
end DO
end if
Do k=1,size(x2,1)
cellk=spldom%cellk(i(k),j(k))
if(cellk%spldirkind.eq.0)then
if(gtmp(k,0) .ge. 0) then
if(size(g,2) .gt. 1) then
g(k,1:2)=-gtmp(k,1:2)*spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val
end if
g(k,0)=(1-gtmp(k,0))*spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val
else
g(k,0)=spldom%boundaries(cellk%linkedboundaries(1))%Dirichlet_val
if(size(g,2).gt. 1) then
g(k,1:2)=0
end if
end if
else
g(k,:)=0
end if
end DO
End SUBROUTINE spline_g
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Evaluates the geometric weight induced by the spline curve defined by b_curve at position (z,r)
!> @param[in] b_curve spline_boundary containing the spline curve parameters
!> @param[in] z axial position where the weight is evaluated
!> @param[in] r radial position where the weight is evaluated
!> @param[out] weight(:) weight index defines the order of derivation by r or z
!> @param[in] h distance from the spline at which the weight is 1
!> @param[out] distance unscaled distance between evaluation point and spline b_curve
!> @param[inout] leftknot initial guess for the closest spline knot of the points (r,z)
!---------------------------------------------------------------------------
subroutine splineweight(b_curve, z, r, weight, h, distance, guess, lguess, rguess)
Use forSISL, ONLY: s1227,s1221, s1774
type(spline_boundary):: b_curve
Real(kind=db)::r,z
Real(kind=db):: weight(0:)
Real(kind=db),OPTIONAL:: distance
real(kind=db),OPTIONAL:: guess
real(kind=db),OPTIONAL:: lguess
real(kind=db),OPTIONAL:: rguess
integer:: sstatus, der, left,siz
real(kind=db):: h, d, tpos, proj
real(kind=real64):: curvepos(2*b_curve%curve%idim)
real(kind=db):: leftpar, rightpar,guesspar
weight=0
der=1
sstatus=-1
guesspar=-1.0_db
if(present(lguess) .and. present(rguess)) then
leftpar=lguess
rightpar=rguess
guesspar=(lguess+rguess)/2
call s1774(b_curve%curve,(/z,r/),b_curve%curve%idim,b_curve%epsge,leftpar,rightpar,guesspar,tpos,sstatus)
if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1774 for splineweight at ", z, r
else
call dist(b_curve,(/z,r/),d,tpos)
end if
! position and derivative wrt r,z
call s1221(b_curve%curve,der,tpos,left,curvepos,sstatus)
if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", z, r
if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", z, r
d=sqrt((curvepos(1)-z)**2+(curvepos(2)-r)**2)
weight(0)=1-max((h-d)/h,0.0_db)**3
curvepos(3:4)=curvepos(3:4)/sqrt(curvepos(3)**2+curvepos(4)**2)
! if the projection of the distance vector on the normal is negative, the weight is negative
proj=(-(z-curvepos(1))*curvepos(4)+(r-curvepos(2))*curvepos(3))
if (proj .lt. 0 .or. abs(abs(proj) -sqrt((z-curvepos(1))**2+(r-curvepos(2))**2)).gt.1e-10) weight(0)=-weight(0)
!if (proj .lt. 0 ) weight(0)=-weight(0)
siz=size(weight,1)
if (size(weight,1).gt.1 .and. abs(weight(0)) .lt. 1) then
weight(1)=-3*curvepos(4)*(h-d)**2/h**3
weight(2)=+3*curvepos(3)*(h-d)**2/h**3
end if
if(present(distance)) distance=d
if(present(guess)) guess=tpos
end subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> Calculates the closest distance between the point and the spline b_curve
!> @param[in] b_curve spline_boundary containing the spline curve parameters
!> @param[in] point(:) array containing the position from which to calculate the distance
!> @param[out] distance distance from the point to the spline
!> @param[in] pos parameter value of the closest point on the spline
!---------------------------------------------------------------------------
subroutine dist(b_curve, point, distance, pos)
Use forSISL, ONLY: s1957,s1953, s1221,s1227
type(spline_boundary):: b_curve
Real(kind=db):: point(:)
real(kind=db):: distance
Real(kind=db),optional::pos
REAL(real64):: posres, epsco, epsge,curvepos(2),d,distmin
REAL(real64),allocatable::intpar(:)
integer:: numintpt, sstat, numintcu,i,left,sstatus
type(SISLIntCurve),ALLOCATABLE:: intcurve(:)
epsco=1.0e-15
epsge=1.0e-15
!epsco=0
!epsge=b_curve%epsge
numintpt=0
sstatus=0
distmin=HUGE(d)
call s1953(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,numintpt,intpar,numintcu,intcurve,sstatus)
if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1953 for splineweight at ", point(1), point(2)
if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1953 for splineweight at ",point(1), point(2)
if(numintpt .gt. 1) then
Do i=1,numintpt
call s1227(b_curve%curve,0,intpar(i),left,curvepos,sstatus)
if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1221 for splineweight at ", point(1), point(2)
if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1221 for splineweight at ",point(1), point(2)
d=(curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2
if(d .lt. distmin) then
distmin=d
posres=intpar(i)
end if
end do
else if(numintpt .gt. 0) then
posres=intpar(1)
end if
if(numintcu.ge.1) then
posres=0.5*(intcurve(1)%epar1(1)+intcurve(1)%epar1(2))
end if
!call s1957(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,posres,distance,sstat)
!if (sstat > 1 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1957 for splineweight at ", point
!if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1957 for splineweight at ", point
if (present(pos)) pos=posres
END subroutine
SUBROUTINE classify(x1, x2, cellk, spldom, wpredir, wpretot)
real(kind=db), INTENT(IN):: x2(2), x1(2)
type(cellkind), intent(INOUT):: cellk
type(spline_domain)::spldom
Real(kind=db):: zeval(4),reval(4), wpretot, wpredir
real(kind=db), allocatable:: guess(:,:), w(:,:,:)
Real(kind=db):: dmin, insidedir, insidetot, distance
integer:: i,k
allocate(guess(spldom%nbsplines,4))
allocate(w(0:2,spldom%nbsplines,4))
w=0
cellk%spldirkind=0
guess=-1.0_db
dmin=HUGE(spldom%dist_extent)
cellk%linkedboundaries=0
zeval=(/ x1(1),x1(2),x1(1),x1(2) /)
reval=(/ x2(1),x2(1),x2(2),x2(2) /)
insidedir=1
do i=1,spldom%nbsplines
do k=1,4
! calculate the weight for each spline boundaries at each cell corner
call splineweight(spldom%boundaries(i),zeval(k),reval(k),w(:,i,k),spldom%dist_extent,distance,guess(i,k))
! We find the closest boundary to this point
if(distance .lt. dmin) then
! If we are close enough we check if we are below dist_extent and need to calculate the distance each time
if(distance .lt. spldom%dist_extent) then
if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then
cellk%linkedboundaries(1)=i
cellk%spldirkind=0
end if
cellk%linkedboundaries(2)=i
cellk%spltotkind=0
end if
dmin=distance
! Otherwise we define the interior by the closest spline
if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then
insidedir=w(0,i,k)
end if
insidetot=w(0,i,k)
end if
end do
end do
if(cellk%linkedboundaries(1) .gt. 0) then
i=cellk%linkedboundaries(1)
cellk%lguess(1)=minval(guess(i,:),1,guess(i,:).ge.0)
cellk%rguess(1)=maxval(guess(i,:),1)
wpredir=w(0,i,1)
else
cellk%spldirkind=sign(1,int(insidedir))
wpredir=insidedir
end if
if(cellk%linkedboundaries(2) .gt. 0) then
i=cellk%linkedboundaries(2)
wpretot=w(0,i,1)
else
cellk%spltotkind=sign(1,int(insidetot))
wpretot=insidetot
end if
end subroutine
subroutine classifycells(spldom)
use forSISL, ONLY: s1537, s1424
type(spline_domain):: spldom
integer:: i,j,sstat
type(cellkind):: cellk
real(kind=db), allocatable:: wpretot(:,:,:), wpredir(:,:,:)
allocate(wpretot(1:1,0:spldom%nb1,0:spldom%nb2))
allocate(wpredir(1:1,0:spldom%nb1,0:spldom%nb2))
wpretot=0
wpredir=0
!$OMP PARALLEL DO private(i,j)
do i=0,spldom%nb1-1
!DIR$ UNROLL
do j=0,spldom%nb2-1
call classify(spldom%x1(i:i+1),spldom%x2(j:j+1),spldom%cellk(i,j),spldom, wpredir(1,i,j),wpretot(1,i,j))
end do
end do
!$OMP END PARALLEL DO
i=spldom%nb1
!$OMP PARALLEL DO private(j,cellk)
do j=0,spldom%nb2
cellk=spldom%cellk(min(i,spldom%nb1-1),min(j,spldom%nb2-1))
If(abs(cellk%spldirkind) .eq. 1) Then
wpredir(1,i,j)=cellk%spldirkind
else
call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), spldom%x1(i),spldom%x2(j), wpredir(:,i,j),spldom%dist_extent)
end IF
If(abs(cellk%spltotkind) .eq. 1) Then
wpretot(1,i,j)=cellk%spltotkind
else
call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), spldom%x1(i),spldom%x2(j), wpretot(:,i,j),spldom%dist_extent)
end IF
end do
!$OMP END PARALLEL DO
!$OMP PARALLEL DO private(i,cellk)
do i=0,spldom%nb1
j=spldom%nb2
cellk=spldom%cellk(min(i,spldom%nb1-1),min(j,spldom%nb2-1))
If(abs(cellk%spldirkind) .eq. 1) Then
wpredir(1,i,j)=cellk%spldirkind
else
call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), spldom%x1(i),spldom%x2(j), wpredir(:,i,j),spldom%dist_extent)
end IF
If(abs(cellk%spltotkind) .eq. 1) Then
wpretot(1,i,j)=cellk%spltotkind
else
call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), spldom%x1(i),spldom%x2(j), wpretot(:,i,j),spldom%dist_extent)
end IF
end do
!$OMP END PARALLEL DO
call s1537(reshape(wpredir(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%Dirdomweight,sstat)
if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight"
if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight"
call s1537(reshape(wpretot(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%totdomweight,sstat)
if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight"
if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight"
end subroutine
subroutine getindex(x1,x2,spldom, i, j)
use distrib, ONLY: closest
type(spline_domain):: spldom
real(kind=db):: x1(:), x2(:)
integer:: i(:),j(:)
call locintv(spldom%splrz%sp1,x1, i)
call locintv(spldom%splrz%sp2,x2, j)
end subroutine
subroutine free_bsplinecurve(b_curve)
type(spline_boundary):: b_curve
call freeCurve(b_curve%curve)
!call freeIntCurve(b_curve%intcurve)
end subroutine
END MODULE splinebound
c4science · Help