diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index e5720b0..8275ef1 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,756 +1,758 @@ 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 + periodic=.false. + 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 forSISL,ONLY: newcurve, s1630 use mpi type(spline_boundary):: b_curve Real(kind=db):: cpoints(:,:) - Real(REAL64),ALLOCATABLE:: knots(:),points(:) + Real(REAL64),ALLOCATABLE:: points(:) + Real(REAL64):: astpar 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 + Integer:: nbpoints,i, dim, jstat, bsptype 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) + !CALL newCurve(nbpoints, order, knots, points, kind, dim, copy, b_curve%curve) + bsptype=1 ! open boundaries b-spline + if(period) bsptype=-1 ! closed periodic curve + astpar=0.0 ! starting parameter for the knots vector + CALL s1630(points, nbpoints, astpar, bsptype, dim, order, b_curve%curve, jstat) + if (jstat > 0 ) WRITE(*,*) "Warning ", jstat," in curve initialisation s1630 for splineweight" + if (jstat < 0 ) WRITE(*,*) "Error ", jstat," in curve initialisation s1630 for splineweight" 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