gradient_reconstruction.F File Reference

#include "implicit_f.inc"
#include "vect01_c.inc"
#include "com04_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	gradient_reconstruction (ixs, x, ale_connect, itrimat, segvar)

Function/Subroutine Documentation

gradient_reconstruction()

subroutine gradient_reconstruction	(	integer, dimension(nixs,numels), intent(in)	ixs,
		dimension(3,numnod), intent(in)	x,
		type(t_ale_connectivity), intent(in)	ale_connect,
		integer, intent(in)	itrimat,
		type(t_segvar)	segvar )

Definition at line 35 of file gradient_reconstruction.F.

C-----------------------------------------------
C  D e s c r i p t i o n
C  This subroutine computes a gradient of the scalar field value in each 
C  element:
C        mean square approximation of the gradient in the face related 
C        neighborhood of the element
C-----------------------------------------------
C   M o d u l e s 
C-----------------------------------------------
      USE alemuscl_mod
      USE segvar_mod
      USE ale_connectivity_mod
      use element_mod , only :nixs
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include "vect01_c.inc"
#include "com04_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: IXS(NIXS,NUMELS)
      my_real, INTENT(IN) :: x(3,numnod)
      INTEGER, INTENT(IN) :: ITRIMAT
      TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I, II, KK, IAD2, LGTH
      my_real :: xk, yk, zk, xl, yl, zl, xf, yf, zf
      my_real :: valk, vall
      my_real :: mat(3, 3), rhs(3), sol(3)
      INTEGER :: VOIS_ID
      INTEGER :: FACE_TO_NODE_LOCAL_ID(6, 4)
      my_real :: norm(3), a(3), b(3), c(3), surf, surf1, surf2
      TYPE(t_segvar) :: SEGVAR
C-----------------------------------------------
C   S o u r c e   L i n e s 
C-----------------------------------------------
!!! Once for all, associate node local id to a face number
!!! Face 1
      face_to_node_local_id(1, 1) = 1 ; face_to_node_local_id(1, 2) = 4
      face_to_node_local_id(1, 3) = 3 ; face_to_node_local_id(1, 4) = 2
!!! Face 2
      face_to_node_local_id(2, 1) = 3 ; face_to_node_local_id(2, 2) = 4
      face_to_node_local_id(2, 3) = 8 ; face_to_node_local_id(2, 4) = 7
!!! Face 3
      face_to_node_local_id(3, 1) = 5 ; face_to_node_local_id(3, 2) = 6
      face_to_node_local_id(3, 3) = 7 ; face_to_node_local_id(3, 4) = 8
!!! Face 4
      face_to_node_local_id(4, 1) = 1 ; face_to_node_local_id(4, 2) = 2
      face_to_node_local_id(4, 3) = 6 ; face_to_node_local_id(4, 4) = 5
!!! Face 5
      face_to_node_local_id(5, 1) = 2 ; face_to_node_local_id(5, 2) = 3
      face_to_node_local_id(5, 3) = 7 ; face_to_node_local_id(5, 4) = 6
!!! Face 6
      face_to_node_local_id(6, 1) = 1 ; face_to_node_local_id(6, 2) = 5
      face_to_node_local_id(6, 3) = 8 ; face_to_node_local_id(6, 4) = 4  
 
      DO i = lft, llt
         ii = i + nft
         !!! Reset mat, rhs
         mat(1:3, 1:3) = zero ; rhs(1:3) = zero
         !!! Value of the target function in the element
         valk = alemuscl_buffer%VOLUME_FRACTION(ii,itrimat)
         xk = alemuscl_buffer%ELCENTER(ii,1) ;
         yk = alemuscl_buffer%ELCENTER(ii,2) ;
         zk = alemuscl_buffer%ELCENTER(ii,3)
         !!! IXS(2:9, II) : Node global ID
         iad2 = ale_connect%ee_connect%iad_connect(ii)
         lgth = ale_connect%ee_connect%iad_connect(ii+1)-iad2
         DO kk = 1, lgth
            vois_id = ale_connect%ee_connect%connected(iad2 + kk - 1)
            !IF (VOIS_ID /= 0 .AND. VOIS_ID <= NUMELS) THEN
            IF (vois_id > 0) THEN
               !!! Value of the target function in the current neighbor
               vall = alemuscl_buffer%VOLUME_FRACTION(vois_id,itrimat)
               xl = alemuscl_buffer%ELCENTER(vois_id,1) ;
               yl = alemuscl_buffer%ELCENTER(vois_id,2) ; 
               zl = alemuscl_buffer%ELCENTER(vois_id,3) ;
            ELSE
               IF(vois_id == 0) THEN
                 vall = valk
               ELSE
                 !vois_id<0 : means EBCS), -vois_id is seg_id
                 vall = segvar%PHASE_ALPHA(itrimat,-vois_id)
               ENDIF
               xf = zero
               yf = zero
               zf = zero
 
               a(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 1) + 1, ii))
               b(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 2) + 1, ii))
               c(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 3) + 1, ii))
 
               norm(1) = (b(2) - a(2)) * (c(3) - a(3)) - (b(3) - a(3)) * (c(2) - a(2))
               norm(2) = (b(3) - a(3)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(3) - a(3))
               norm(3) = (b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1))
 
               surf1 = half * abs(sqrt(norm(1) * norm(1) + norm(2) * norm(2) + norm(3) * norm(3)))
               xf = surf1 * third * (a(1) + b(1) + c(1))
               yf = surf1 * third * (a(2) + b(2) + c(2))
               zf = surf1 * third * (a(3) + b(3) + c(3))
 
               a(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 1) + 1, ii))
               b(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 3) + 1, ii))
               c(1:3) = x(1:3, ixs(face_to_node_local_id(kk, 4) + 1, ii))
 
               norm(1) = (b(2) - a(2)) * (c(3) - a(3)) - (b(3) - a(3)) * (c(2) - a(2))
               norm(2) = (b(3) - a(3)) * (c(1) - a(1)) - (b(1) - a(1)) * (c(3) - a(3))
               norm(3) = (b(1) - a(1)) * (c(2) - a(2)) - (b(2) - a(2)) * (c(1) - a(1))
 
               surf2 = half * abs(sqrt(norm(1) * norm(1) + norm(2) * norm(2) + norm(3) * norm(3)))
               xf = xf + surf2 * third * (a(1) + b(1) + c(1))
               yf = yf + surf2 * third * (a(2) + b(2) + c(2))
               zf = zf + surf2 * third * (a(3) + b(3) + c(3)) 
               
               surf = surf1 + surf2
               xf = xf / surf
               yf = yf / surf
               zf = zf / surf
               !!! Build face centroid
c$$$               DO JJ = 1, 4
c$$$                  NODE_ID = IXS(FACE_TO_NODE_LOCAL_ID(KK, JJ) + 1, II)
c$$$                  XF = XF + FOURTH * X(1, NODE_ID)
c$$$                  YF = YF + FOURTH * X(2, NODE_ID)
c$$$                  ZF = ZF + FOURTH * X(3, NODE_ID)
c$$$               ENDDO
               xl = two * xf - alemuscl_buffer%ELCENTER(ii,1)
               yl = two * yf - alemuscl_buffer%ELCENTER(ii,2)
               zl = two * zf - alemuscl_buffer%ELCENTER(ii,3)
            ENDIF
 
            !!! Incrementing mat and rhs
            rhs(1) = rhs(1) + (valk - vall) * (xl - xk)
            rhs(2) = rhs(2) + (valk - vall) * (yl - yk)
            rhs(3) = rhs(3) + (valk - vall) * (zl - zk)
            mat(1, 1) = mat(1, 1) + (xl - xk) * (xl - xk)
            mat(1, 2) = mat(1, 2) + (xl - xk) * (yl - yk)
            mat(1, 3) = mat(1, 3) + (xl - xk) * (zl - zk)
            mat(2, 1) = mat(2, 1) + (yl - yk) * (xl - xk)
            mat(2, 2) = mat(2, 2) + (yl - yk) * (yl - yk)
            mat(2, 3) = mat(2, 3) + (yl - yk) * (zl - zk)
            mat(3, 1) = mat(3, 1) + (zl - zk) * (xl - xk)
            mat(3, 2) = mat(3, 2) + (zl - zk) * (yl - yk)
            mat(3, 3) = mat(3, 3) + (zl - zk) * (zl - zk)
         ENDDO
 
         CALL cg(3, mat, rhs, sol, 3, em10)
         !CALL DIRECT_SOLVE(mat, rhs ,sol)
         !!! Solution goes to the gradient
         alemuscl_buffer%GRAD(ii,1,itrimat) = -sol(1)
         alemuscl_buffer%GRAD(ii,2,itrimat) = -sol(2)
         alemuscl_buffer%GRAD(ii,3,itrimat) = -sol(3)
      ENDDO  ! I = LFT, LLT
C-----------------------------------------------