interpolateinface.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	interpolateinface (kk, xstate, xstateini, numpts, nstate_, mi, islavsurf, pslavsurf, ne0, islavsurfold, pslavsurfold)

Function/Subroutine Documentation

interpolateinface()

subroutine interpolateinface	(	integer	kk,
		real*8, dimension(nstate_,mi(1),*)	xstate,
		real*8, dimension(nstate_,mi(1),*)	xstateini,
		integer	numpts,
		integer	nstate_,
		integer, dimension(*)	mi,
		integer, dimension(2,*)	islavsurf,
		real*8, dimension(3,*)	pslavsurf,
		integer	ne0,
		integer, dimension(2,*)	islavsurfold,
		real*8, dimension(3,*)	pslavsurfold
	)

!
      implicit none
!
      integer numpts,i_int,n_int,i,j,k,kk,l,ll,nn,
     &  koncont(3,2*numpts+1),itri,kflag,neigh(1),kneigh,
     &  imastop(3,2*numpts+1),indexcj,nopespringj,list(numpts),
     &  igauss,mi(*),nstate_,itriangle(100),itriold,
     &  ifaceq(8,6),ip(numpts),ne0,iloc,itrinew,ntriangle,
     &  ifacet(6,4),ifacew1(4,5),ifacew2(8,5),n,islavsurf(2,*),
     &  ibin(numpts),ivert1,ntriangle_,nterms,m,islavsurfold(2,*),
     &  nx(2*numpts+1),ny(2*numpts+1),isol,id
!
      real*8 xstate(nstate_,mi(1),*),p(3),pslavsurfold(3,*),
     &  xstateini(nstate_,mi(1),*),coi(2,numpts+3),pneigh(3,3),
     &  cg(2,2*numpts+1),xdist,pslavsurf(3,*),xil,etl,ratio(3),
     &  z(9,3),x(2*numpts+1),xo(2*numpts+1),
     &  y(2*numpts+1),yo(2*numpts+1),straight(9,2*numpts+1),dist
!
!     nodes per face for hex elements
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
!
!     nodes per face for tet elements
!
      data ifacet /1,3,2,7,6,5,
     &             1,2,4,5,9,8,
     &             2,3,4,6,10,9,
     &             1,4,3,8,10,7/
!
!     nodes per face for linear wedge elements
!
      data ifacew1 /1,3,2,0,
     &             4,5,6,0,
     &             1,2,5,4,
     &             2,3,6,5,
     &             3,1,4,6/
!
!     nodes per face for quadratic wedge elements
!
      data ifacew2 /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             3,1,4,6,9,13,12,15/
!
      kneigh=1
!
      do i=1,numpts
         list(i)=i
      enddo
!
!       Loop over the old integration points within the face
!
      ll=0
      do l=islavsurfold(2,kk)+1,islavsurfold(2,kk+1)
         ll=ll+1
         coi(1,ll)=pslavsurfold(1,l)
         coi(2,ll)=pslavsurfold(2,l)
         ip(ll)=ne0+l
      enddo
!
!     Calling deltri for triangulation
!     
      kflag=2
!     
      call deltri(numpts,numpts,coi(1,1:numpts+3),coi(2,1:numpts+3),
     &            list,ibin,koncont,imastop,n)
!
!       Arranging imastop field corresponding to "imastop"
!
      nn=0
      do i=1,n
         ivert1=imastop(1,i)
         imastop(1,i)=imastop(2,i)
         imastop(2,i)=imastop(3,i)
         imastop(3,i)=ivert1
      enddo
!
!     determining the center of gravity and the bounding planes
!     of the triangles
!
      call updatecont2d(koncont,n,coi,cg,straight)
!    
!     sorting the centers of gravity
! 
      do l=1,n
         xo(l)=cg(1,l)
         x(l)=xo(l)
         nx(l)=l
         yo(l)=cg(2,l)
         y(l)=yo(l)
         ny(l)=l
      enddo
      kflag=2
      call dsort(x,nx,n,kflag)
      call dsort(y,ny,n,kflag)
!
!       Loop over the new integration points
!
      do iloc=islavsurf(2,kk)+1,islavsurf(2,kk+1)
         xdist=9.d0
         igauss=iloc
!
!        coordinates of the new integration point
!
         p(1)=pslavsurf(1,igauss)
         p(2)=pslavsurf(2,igauss)
!
!        closest triangle center of gravity
!     
         call near2d(xo,yo,x,y,nx,ny,p(1),p(2),
     &        n,neigh,kneigh)
!     
         isol=0
!     
         itriold=0
         itri=neigh(1)
         ntriangle=0
         ntriangle_=100
!     
         loop1: do
            do l=1,3
               ll=3*l-2
               dist=straight(ll,itri)*p(1)+
     &              straight(ll+1,itri)*p(2)+
     &              straight(ll+2,itri)
!     
               if(dist.gt.0.d0) then
                  itrinew=imastop(l,itri)
                  if(itrinew.eq.0) then
c     write(*,*) '**border reached'
                     exit loop1
                  elseif(itrinew.eq.itriold) then
c     write(*,*) '**solution in between triangles'
                     isol=itri
                     exit loop1
                  else
                     call nident(itriangle,itrinew,ntriangle,id)
                     if(id.gt.0) then
                        if(itriangle(id).eq.itrinew) then
c     write(*,*) '**circular path;no solution'
                           exit loop1
                        endif
                     endif
                     ntriangle=ntriangle+1
                     if(ntriangle.gt.ntriangle_) then
c     write(*,*) '**too many iterations'
                        exit loop1
                     endif
                     do k=ntriangle,id+2,-1
                        itriangle(k)=itriangle(k-1)
                     enddo
                     itriangle(id+1)=itrinew
                     itriold=itri
                     itri=itrinew
                     cycle loop1
                  endif
               elseif(l.eq.3) then
c     write(*,*) '**regular solution'
                  isol=itri
                  exit loop1
               endif
            enddo
         enddo loop1
!
         if(isol.eq.0) then
!
!           mapping the new integration point onto the nearest 
!           triangle
!
            do k=1,3
               do m=1,2
                  pneigh(m,k)=coi(m,koncont(k,itri))
               enddo
               pneigh(3,k)=0.d0
            enddo
            p(3)=0.d0
            nterms=3
!
            call attach(pneigh,p,nterms,ratio,dist,xil,etl)
!
            do m=1,2
               p(m)=0.d0
               do k=1,3
                  p(m)=p(m)+ratio(k)*pneigh(m,k)
               enddo
            enddo
         endif
!     
!     Assigning xstate values from the vertices of the triangle
!     to the field z (integration point values from xstateini)
!     
         do k=1,3
            do i=1,9
               z(i,k)=xstateini(i,1,ip(koncont(k,itri)))
            enddo
         enddo
!     
!     Calling plane_eq to interpolate values using plane equation
!     
         do i=1,9   
            call plane_eq(
     &           coi(1,koncont(1,itri)),coi(2,koncont(1,itri)),z(i,1),
     &           coi(1,koncont(2,itri)),coi(2,koncont(2,itri)),z(i,2),
     &           coi(1,koncont(3,itri)),coi(2,koncont(3,itri)),z(i,3),
     &           p(1),p(2),xstate(i,1,ne0+iloc))
         enddo
!     
      enddo
!     
      return