NEWTON_RAPSON.f90 File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	newton_rapson (xnod, ynod, znod, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, tt, csi_s, eta_s, zeta_s, nofi, id_mpi, ipiu, imeno, toll1, toll2, flag)
	Performs the Newton Rapson algorithm.

Function/Subroutine Documentation

newton_rapson()

subroutine newton_rapson	(	real*8, intent(in)	xnod,
		real*8, intent(in)	ynod,
		real*8, intent(in)	znod,
		real*8	alfa11,
		real*8	alfa12,
		real*8	alfa13,
		real*8	alfa21,
		real*8	alfa22,
		real*8	alfa23,
		real*8	alfa31,
		real*8	alfa32,
		real*8	alfa33,
		real*8	beta11,
		real*8	beta12,
		real*8	beta13,
		real*8	beta21,
		real*8	beta22,
		real*8	beta23,
		real*8	beta31,
		real*8	beta32,
		real*8	beta33,
		real*8	gamma1,
		real*8	gamma2,
		real*8	gamma3,
		real*8	delta1,
		real*8	delta2,
		real*8	delta3,
		integer*4, intent(out)	tt,
		real*8, intent(out)	csi_s,
		real*8, intent(out)	eta_s,
		real*8, intent(out)	zeta_s,
		integer*4, intent(in)	nofi,
		integer*4, intent(in)	id_mpi,
		integer*4, intent(in)	ipiu,
		integer*4, intent(in)	imeno,
		real*8, intent(in)	toll1,
		real*8, intent(in)	toll2,
		integer*4, intent(in)	flag
	)

Performs the Newton Rapson algorithm.

Author

Ilario Mazzieri

Date

September, 2013

Version

1.0

Parameters

[in]	xnod	x-coord of the point
[in]	ynod	y-coord of the point
[in]	znod	z-coord of the point
[in]	alfa11	costant for the bilinear map
[in]	alfa12	costant for the bilinear map
[in]	alfa13	costant for the bilinear map
[in]	alfa21	costant for the bilinear map
[in]	alfa22	costant for the bilinear map
[in]	alfa23	costant for the bilinear map
[in]	alfa31	costant for the bilinear map
[in]	alfa32	costant for the bilinear map
[in]	alfa33	costant for the bilinear map
[in]	beta11	costant for the bilinear map
[in]	beta12	costant for the bilinear map
[in]	beta13	costant for the bilinear map
[in]	beta21	costant for the bilinear map
[in]	beta22	costant for the bilinear map
[in]	beta23	costant for the bilinear map
[in]	beta31	costant for the bilinear map
[in]	beta32	costant for the bilinear map
[in]	beta33	costant for the bilinear map
[in]	gamma1	costant for the bilinear map
[in]	gamma2	costant for the bilinear map
[in]	gamma3	costant for the bilinear map
[in]	delta1	costant for the bilinear map
[in]	delta2	costant for the bilinear map
[in]	delta3	costant for the bilinear map
[in]	tt	1 if the node belongs to the element, 0 otherwise
[in]	nofi	number of iterations (dummy)
[in]	id_mpi	mpi processor identity
[in]	ipiu	control parameter (dummy)
[in]	imeno	control parameter (dummy)
[in]	toll1	fixed tolerance
[in]	toll2	fixed tolerance
[in]	flag	control parameter
[out]	csi_s	x-coordinate in the reference element of xnod,ynod,znod
[out]	eta_s	y-coordinate in the reference element of xnod,ynod,znod
[out]	zeta_s	z-coordinate in the reference element of xnod,ynod,znod

Definition at line 62 of file NEWTON_RAPSON.f90.

 
      implicit none
                      
      integer*4 :: nofiter
      integer*4 :: per_csi, per_eta, per_zeta
      integer*4 ::  iic, il,ih
      integer*4, intent(out) :: tt 
      integer*4, intent(in) :: nofi,id_mpi,ipiu,imeno, flag
      
      real*8 :: ll,lm ,xvt1,yvt1,xvt2,yvt2,xvc1,yvc1,xvc2,yvc2
      real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3
      real*8 :: e1,e2,e3,f1,f2,f3,g1,g2,g3,h1,h2,h3
      real*8 :: delta_csi, delta_eta, delta_zeta
      real*8 :: a,b,c,d,e,f,g,h,p, det
      real*8 :: a11, a12, a13, a21, a22, a23, a31, a32, a33
      real*8 :: ff1, ff2, ff3
      real*8 :: alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,alfa31,alfa32,alfa33
      real*8 :: beta11,beta12,beta13,beta21,beta22,beta23,beta31,beta32,beta33
      real*8 :: gamma1,gamma2,gamma3,delta1,delta2,delta3
      real*8, intent(in) :: xnod, ynod, znod, toll1, toll2
      real*8, intent(out) ::  csi_s, eta_s, zeta_s
 
 
      a1 = 8.d0*gamma1
      a2 = 8.d0*gamma2
      a3 = 8.d0*gamma3
      b1 = 8.d0*beta13
      b2 = 8.d0*beta23
      b3 = 8.d0*beta33
      c1 = 8.d0*beta11
      c2 = 8.d0*beta21
      c3 = 8.d0*beta31
      d1 = 8.d0*beta12
      d2 = 8.d0*beta22
      d3 = 8.d0*beta32
      e1 = 8.d0*alfa11
      e2 = 8.d0*alfa21
      e3 = 8.d0*alfa31
      f1 = 8.d0*alfa12
      f2 = 8.d0*alfa22
      f3 = 8.d0*alfa32
      g1 = 8.d0*alfa13
      g2 = 8.d0*alfa23
      g3 = 8.d0*alfa33
      h1 = 8.d0*delta1
      h2 = 8.d0*delta2
      h3 = 8.d0*delta3
 
      csi_s = 0.d0
      eta_s = 0.d0
      zeta_s = 0.d0
      delta_csi = 0.0d0
      delta_eta = 0.0d0 
      delta_zeta = 0.0d0     
      
      tt = 0
      nofiter = 2000
 
      per_csi = 0
      per_eta = 0
      per_zeta = 0
 
 
      do il = 1, nofiter
                
        a = a1*eta_s*zeta_s + b1*eta_s + d1*zeta_s + e1
        b = a1*csi_s*zeta_s + b1*csi_s + c1*zeta_s + f1
        c = a1*csi_s*eta_s + c1*eta_s + d1*csi_s + g1
 
        d = a2*eta_s*zeta_s + b2*eta_s + d2*zeta_s + e2
        e = a2*csi_s*zeta_s + b2*csi_s + c2*zeta_s + f2
        f = a2*csi_s*eta_s + c2*eta_s + d2*csi_s + g2
 
        g = a3*eta_s*zeta_s + b3*eta_s + d3*zeta_s + e3
        h = a3*csi_s*zeta_s + b3*csi_s + c3*zeta_s + f3 
        p = a3*csi_s*eta_s + c3*eta_s + d3*csi_s + g3
 
        det = a*(e*p - f*h) - b*(d*p - f*g) + c*(d*h -e*g)
        
        if(det .eq. 0.d0) then
          return
        endif
 
        a11 = e*p - f*h
        a12 = -d*p + f*g
        a13 = d*h - e*g
        a21 = -b*p + c*h 
        a22 = a*p - c*g
        a23 = -a*h + b*g
        a31 = b*f - c*e
        a32 = -a*f + c*d
        a33 = a*e - b*d
        
        ff1 = a1*csi_s*eta_s*zeta_s + b1*csi_s*eta_s + c1*eta_s*zeta_s + d1*csi_s*zeta_s &
                                   + e1*csi_s + f1*eta_s + g1*zeta_s + h1 -8.d0*xnod 
        
        ff2 = a2*csi_s*eta_s*zeta_s + b2*csi_s*eta_s + c2*eta_s*zeta_s + d2*csi_s*zeta_s &
                                   + e2*csi_s + f2*eta_s + g2*zeta_s + h2 -8.d0*ynod 
        
        ff3 = a3*csi_s*eta_s*zeta_s + b3*csi_s*eta_s + c3*eta_s*zeta_s + d3*csi_s*zeta_s &
                                   + e3*csi_s + f3*eta_s + g3*zeta_s + h3 -8.d0*znod 
 
 
        delta_csi = (-1.d0/det)*(a11*ff1 + a21*ff2 + a31*ff3)
 
        delta_eta = (-1.d0/det)*(a12*ff1 + a22*ff2 + a32*ff3)
 
        delta_zeta = (-1.d0/det)*(a13*ff1 + a23*ff2 + a33*ff3)
 
 
         
        if (abs(delta_csi).le. toll1 .and. abs(delta_eta).le. toll1 .and. abs(delta_zeta).le. toll1 ) then
            tt = 1
           
           if(flag.eq.1) then 
              if (dabs(csi_s) .gt. toll2 .or. dabs(eta_s) .gt. toll2 .or. dabs(zeta_s).gt. toll2 ) then
              tt = 0
              return
            endif
          endif
         
          return
        endif
        
       
        csi_s = csi_s + delta_csi
        eta_s = eta_s + delta_eta
        zeta_s = zeta_s + delta_zeta
        
        if (csi_s .gt. 1.d0 .AND. csi_s .lt. 2.d0 .AND. per_csi .eq. 0) then
            csi_s = 1.d0 
            per_csi = 1
        endif
             
        if (csi_s .lt. -1.d0 .AND. csi_s .gt. -2.d0 .AND. per_csi .eq. 0) then
            csi_s = -1.d0
            per_csi = 1
        endif
        
        
        if (eta_s .gt. 1.d0 .AND. eta_s .lt. 2.d0 .AND. per_eta .eq. 0) then
            eta_s = 1.d0 
            per_eta = 1
        endif
            
        if (eta_s .lt. -1.d0 .AND. eta_s .gt. -2.d0 .AND. per_eta .eq. 0) then
            eta_s = -1.d0
            per_eta = 1
        endif    
            
        if (zeta_s .gt. 1.d0 .AND. zeta_s .lt. 2.d0 .AND. per_zeta .eq. 0) then
            zeta_s = 1.d0 
            per_zeta = 1
        endif
            
        if (zeta_s .lt. -1.d0 .AND. zeta_s .gt. -2.d0 .AND. per_zeta .eq. 0) then
            zeta_s = -1.d0
            per_zeta = 1
        endif    
 
 
      enddo   
      
      if(il .ge. nofiter) then      
        tt = 0
        return
      endif
 

Referenced by make_dg_interface(), and write_file_dgfs().

Here is the caller graph for this function: