utilities.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_utilities
  use hecmw
  implicit none
 
  real(kind=kreal), parameter, private :: pi=3.14159265358979d0
 
contains
 
  subroutine memget(var,dimn,syze)
    integer :: var,dimn,syze,bite
    parameter(bite=1)
    var = var + dimn*syze*bite
  end subroutine memget
 
 
  subroutine append_int2name( n, fname, n1 )
    integer, intent(in)             :: n
    integer, intent(in), optional  ::  n1
    character(len=*), intent(inout) :: fname
    integer            :: npos, nlen
    character(len=128) :: tmpname, tmp
 
    npos = scan( fname, '.')
    nlen = len_trim( fname )
    if( nlen>128 ) stop "String too long(>128) in append_int2name"
    if( n>100000 ) stop "Integer too big>100000 in append_int2name"
    tmpname = fname
    if( npos==0 ) then
      write( fname, '(a,i6)') fname(1:nlen),n
    else
      write( tmp, '(i6,a)') n,tmpname(npos:nlen)
      fname = tmpname(1:npos-1) // adjustl(tmp)
    endif
    if(present(n1).and.n1/=0)then
      write(tmp,'(i8)')n1
      fname = fname(1:len_trim(fname))//'.'//adjustl(tmp)
    endif
  end subroutine
 
  subroutine insert_int2array( iin, carray )
    integer, intent(in) :: iin
    integer, pointer :: carray(:)
 
    integer :: i, oldsize
    integer, pointer :: dumarray(:) => null()
    if( .not. associated(carray) ) then
      allocate( carray(1) )
      carray(1) = iin
    else
      oldsize = size( carray )
      allocate( dumarray(oldsize) )
      do i=1,oldsize
        dumarray(i) = carray(i)
      enddo
      deallocate( carray )
      allocate( carray(oldsize+1) )
      do i=1,oldsize
        carray(i) = dumarray(i)
      enddo
      carray(oldsize+1) = iin
    endif
    if( associated(dumarray) ) deallocate( dumarray )
  end subroutine
 
  subroutine tensor_eigen3( tensor, eigval )
    real(kind=kreal), intent(in)  :: tensor(6)          
    real(kind=kreal), intent(out) :: eigval(3)     
 
    real(kind=kreal) :: i1,i2,i3,r,theta,q, x(3,3), xx(3,3)
    real(kind=kreal), parameter :: small = 1.0d-10
 
    x(1,1)=tensor(1); x(2,2)=tensor(2); x(3,3)=tensor(3)
    x(1,2)=tensor(4); x(2,1)=x(1,2)
    x(2,3)=tensor(5); x(3,2)=x(2,3)
    x(3,1)=tensor(6); x(1,3)=x(3,1)
 
    xx= matmul( x,x )
    i1= x(1,1)+x(2,2)+x(3,3)
    i2= 0.5d0*( i1*i1 - (xx(1,1)+xx(2,2)+xx(3,3)) )
    i3= x(1,1)*x(2,2)*x(3,3)+x(2,1)*x(3,2)*x(1,3)+x(3,1)*x(1,2)*x(2,3)    &
      -x(3,1)*x(2,2)*x(1,3)-x(2,1)*x(1,2)*x(3,3)-x(1,1)*x(3,2)*x(2,3)
 
    r=(-2.d0*i1*i1*i1+9.d0*i1*i2-27.d0*i3)/54.d0
    q=(i1*i1-3.d0*i2)/9.d0
 
    ! Check discriminant to ensure all eigenvalues are real
    ! For a 3x3 symmetric matrix, Delta = R^2 - Q^3 must be negative or zero
    if(r*r - q*q*q > small) then
      write(*,*) 'ERROR in tensor_eigen3: Discriminant R^2 - Q^3 > 0'
      write(*,*) 'This indicates the algorithm has numerical issues.'
      write(*,*) 'Discriminant =', r*r - q*q*q
      stop 'Eigenvalue calculation failed - possible numerical error'
    endif
 
    ! Add handling for special cases (when Q is close to 0)
    if(abs(q) < small) then
      ! Matrix has nearly equal diagonal elements (e.g. close to identity matrix)
      eigval(1) = i1/3.d0
      eigval(2) = i1/3.d0
      eigval(3) = i1/3.d0
    else
      ! For Q^3 sufficiently small, address numerical stability
      if(abs(q*q*q) < small) then
        ! Q^3 is very small, but not zero - special handling
        if(abs(r) < small) then
          ! Both R and Q^3 small means nearly triple root
          eigval(1) = i1/3.d0
          eigval(2) = i1/3.d0
          eigval(3) = i1/3.d0
        else
          ! When Q^3 is very small, determine theta based on sign of R only
          ! This avoids division by a very small number which could cause instability
          if(r >= 0.0d0) then
            theta = 0.0d0      ! R > 0 case
          else
            theta = pi         ! R < 0 case
          endif
        endif
      else
        ! Ensure R/sqrt(Q^3) is within [-1,1] for numerical stability
        theta = acos(max(-1.0d0, min(1.0d0, r/dsqrt(q*q*q))))
      endif
 
      eigval(1) = -2.d0*dsqrt(q)*cos(theta/3.d0)+i1/3.d0
      eigval(2) = -2.d0*dsqrt(q)*cos((theta+2.d0*pi)/3.d0)+i1/3.d0
      eigval(3) = -2.d0*dsqrt(q)*cos((theta-2.d0*pi)/3.d0)+i1/3.d0
    endif
 
  end subroutine
 
  subroutine eigen3 (tensor, eigval, princ)
    real(kind=kreal), intent(in)  :: tensor(6)     
    real(kind=kreal), intent(out) :: eigval(3)     
    real(kind=kreal), intent(out) :: princ(3, 3)   
 
    integer, parameter :: msweep = 50
    integer :: i,j, is, ip, iq, ir
    real(kind=kreal) :: fsum, od, theta, t, c, s, tau, g, h, hd, btens(3,3), factor
 
    btens(1,1)=tensor(1); btens(2,2)=tensor(2); btens(3,3)=tensor(3)
    btens(1,2)=tensor(4); btens(2,1)=btens(1,2)
    btens(2,3)=tensor(5); btens(3,2)=btens(2,3)
    btens(3,1)=tensor(6); btens(1,3)=btens(3,1)
    !
    !     Initialise princ to the identity
    !
    factor = 0.d0
    do i = 1, 3
      do j = 1, 3
        princ(i, j) = 0.d0
      end do
      princ(i, i) = 1.d0
      eigval(i) = btens(i, i)
      factor = factor + dabs(btens(i,i))
    end do
    !     Scaling and iszero/isnan exception
    if( factor == 0.d0 .or. factor /= factor ) then
      return
    else
      eigval(1:3) = eigval(1:3)/factor
      btens(1:3,1:3) = btens(1:3,1:3)/factor
    end if
 
    !
    !     Starts sweeping.
    !
    do is = 1, msweep
      fsum = 0.d0
      do ip = 1, 2
        do iq = ip + 1, 3
          fsum = fsum + abs( btens(ip, iq) )
        end do
      end do
      !
      !     If the fsum of off-diagonal terms is zero returns
      !
      if ( fsum < 1.d-10 ) then
        eigval(1:3) = eigval(1:3)*factor
        return
      endif
      !
      !     Performs the sweep in three rotations. One per off diagonal term
      !
      do ip = 1, 2
        do iq = ip + 1, 3
          od = 100.d0 * abs(btens(ip, iq) )
          if ( (od+abs(eigval(ip) ) /= abs(eigval(ip) )) &
              .and. (od+abs(eigval(iq) ) /= abs(eigval(iq) ))) then
            hd = eigval(iq) - eigval(ip)
            !
            !    Evaluates the rotation angle
            !
            if ( abs(hd) + od == abs(hd)  ) then
              t = btens(ip, iq) / hd
            else
              theta = 0.5d0 * hd / btens(ip, iq)
              t = 1.d0 / (abs(theta) + sqrt(1.d0 + theta**2) )
              if ( theta < 0.d0 ) t = - t
            end if
            !
            !     Re-evaluates the diagonal terms
            !
            c = 1.d0 / sqrt(1.d0 + t**2)
            s = t * c
            tau = s / (1.d0 + c)
            h = t * btens(ip, iq)
            eigval(ip) = eigval(ip) - h
            eigval(iq) = eigval(iq) + h
            !
            !     Re-evaluates the remaining off-diagonal terms
            !
            ir = 6 - ip - iq
            g = btens(min(ir, ip), max(ir, ip) )
            h = btens(min(ir, iq), max(ir, iq) )
            btens(min(ir, ip), max(ir, ip) ) = g &
              - s * (h + g * tau)
            btens(min(ir, iq), max(ir, iq) ) = h &
              + s * (g - h * tau)
            !
            !     Rotates the eigenvectors
            !
            do ir = 1, 3
              g = princ(ir, ip)
              h = princ(ir, iq)
              princ(ir, ip) = g - s * (h + g * tau)
              princ(ir, iq) = h + s * (g - h * tau)
            end do
          end if
          btens(ip, iq) = 0.d0
        end do
      end do
    end do ! over sweeps
    !
    !     If convergence is not achieved stops
    !
    stop       ' Jacobi iteration unable to converge'
  end subroutine eigen3
 
  real(kind=kreal) function determinant( mat )
    real(kind=kreal) :: mat(6)     
    real(kind=kreal) :: xj(3,3)
 
    xj(1,1)=mat(1); xj(2,2)=mat(2); xj(3,3)=mat(3)
    xj(1,2)=mat(4); xj(2,1)=xj(1,2)
    xj(2,3)=mat(5); xj(3,2)=xj(2,3)
    xj(3,1)=mat(6); xj(1,3)=xj(3,1)
 
    determinant=xj(1,1)*xj(2,2)*xj(3,3)               &
      +xj(2,1)*xj(3,2)*xj(1,3)                   &
      +xj(3,1)*xj(1,2)*xj(2,3)                   &
      -xj(3,1)*xj(2,2)*xj(1,3)                   &
      -xj(2,1)*xj(1,2)*xj(3,3)                   &
      -xj(1,1)*xj(3,2)*xj(2,3)
  end function determinant
 
  real(kind=kreal) function determinant33( XJ )
    real(kind=kreal) :: xj(3,3)
 
    determinant33=xj(1,1)*xj(2,2)*xj(3,3)               &
      +xj(2,1)*xj(3,2)*xj(1,3)                   &
      +xj(3,1)*xj(1,2)*xj(2,3)                   &
      -xj(3,1)*xj(2,2)*xj(1,3)                   &
      -xj(2,1)*xj(1,2)*xj(3,3)                   &
      -xj(1,1)*xj(3,2)*xj(2,3)
  end function determinant33
 
  subroutine fstr_chk_alloc( imsg, sub_name, ierr )
    use hecmw
    character(*) :: sub_name
    integer(kind=kint) :: imsg
    integer(kind=kint) :: ierr
 
    if( ierr /= 0 ) then
      write(imsg,*) 'Memory overflow at ', sub_name
      write(*,*) 'Memory overflow at ', sub_name
      call hecmw_abort( hecmw_comm_get_comm( ) )
    endif
  end subroutine fstr_chk_alloc
 
  subroutine calinverse(NN, A)
    integer, intent(in)             :: NN
    real(kind=kreal), intent(inout) :: a(nn,nn)
 
    integer          :: I, J,K,IW,LR,IP(NN)
    real(kind=kreal) :: w,wmax,pivot,api,eps,det
    data eps/1.0e-35/
    det=1.d0
    lr = 0
    do i=1,nn
      ip(i)=i
    enddo
    do k=1,nn
      wmax=0.d0
      do i=k,nn
        w=dabs(a(i,k))
        if (w.GT.wmax) then
          wmax=w
          lr=i
        endif
      enddo
      pivot=a(lr,k)
      api=abs(pivot)
      if(api.LE.eps) then
        write(*,'(''PIVOT ERROR AT'',I5)') k
        stop
      end if
      det=det*pivot
      if (lr.NE.k) then
        det=-det
        iw=ip(k)
        ip(k)=ip(lr)
        ip(lr)=iw
        do j=1,nn
          w=a(k,j)
          a(k,j)=a(lr,j)
          a(lr,j)=w
        enddo
      endif
      do i=1,nn
        a(k,i)=a(k,i)/pivot
      enddo
      do i=1,nn
        if (i.NE.k) then
          w=a(i,k)
          if (w.NE.0.) then
            do j=1,nn
              if (j.NE.k) a(i,j)=a(i,j)-w*a(k,j)
            enddo
            a(i,k)=-w/pivot
          endif
        endif
      enddo
      a(k,k)=1.d0/pivot
    enddo
 
    do i=1,nn
      k=ip(i)
      if (k.NE.i) then
        iw=ip(k)
        ip(k)=ip(i)
        ip(i)=iw
        do j=1,nn
          w=a(j,i)
          a(j,i)=a(j,k)
          a(j,k)=w
        enddo
      endif
    enddo
 
  end subroutine calinverse
 
  subroutine calsolve(NN, A, b)
    integer, intent(in)             :: NN
    real(kind=kreal), intent(inout) :: a(nn,nn)
    real(kind=kreal), intent(inout) :: b(nn)
 
    integer          :: i, j, k, imax
    real(kind=kreal) :: w, piv
 
    ! Forward elimination
    do k = 1, nn-1
      imax = k
      do i = k+1, nn
        if( dabs(a(i,k)) > dabs(a(imax,k)) ) imax = i
      end do
      if( imax /= k ) then
        do j = k, nn
          w = a(k,j);  a(k,j) = a(imax,j);  a(imax,j) = w
        end do
        w = b(k);  b(k) = b(imax);  b(imax) = w
      end if
      piv = a(k,k)
      do i = k+1, nn
        w = a(i,k) / piv
        do j = k+1, nn
          a(i,j) = a(i,j) - w * a(k,j)
        end do
        b(i) = b(i) - w * b(k)
      end do
    end do
    ! Back substitution
    do i = nn, 1, -1
      do j = i+1, nn
        b(i) = b(i) - a(i,j) * b(j)
      end do
      b(i) = b(i) / a(i,i)
    end do
  end subroutine calsolve
 
  subroutine cross_product(v1,v2,vn)
    real(kind=kreal),intent(in)  ::  v1(3),v2(3)
    real(kind=kreal),intent(out)  ::  vn(3)
 
    vn(1) = v1(2)*v2(3) - v1(3)*v2(2)
    vn(2) = v1(3)*v2(1) - v1(1)*v2(3)
    vn(3) = v1(1)*v2(2) - v1(2)*v2(1)
  end subroutine cross_product
 
  subroutine transformation(jacob, tm)
    real(kind=kreal),intent(in)  ::  jacob(3,3)   
    real(kind=kreal),intent(out)  ::  tm(6,6)      
 
    integer    ::  i,j
 
    do i=1,3
      do j=1,3
        tm(i,j)= jacob(i,j)*jacob(i,j)
      enddo
      tm(i,4) = jacob(i,1)*jacob(i,2)
      tm(i,5) = jacob(i,2)*jacob(i,3)
      tm(i,6) = jacob(i,3)*jacob(i,1)
    enddo
    tm(4,1) = 2.d0*jacob(1,1)*jacob(2,1)
    tm(5,1) = 2.d0*jacob(2,1)*jacob(3,1)
    tm(6,1) = 2.d0*jacob(3,1)*jacob(1,1)
    tm(4,2) = 2.d0*jacob(1,2)*jacob(2,2)
    tm(5,2) = 2.d0*jacob(2,2)*jacob(3,2)
    tm(6,2) = 2.d0*jacob(3,2)*jacob(1,2)
    tm(4,3) = 2.d0*jacob(1,3)*jacob(2,3)
    tm(5,3) = 2.d0*jacob(2,3)*jacob(3,3)
    tm(6,3) = 2.d0*jacob(3,3)*jacob(1,3)
    tm(4,4) = jacob(1,1)*jacob(2,2) + jacob(1,2)*jacob(2,1)
    tm(5,4) = jacob(2,1)*jacob(3,2) + jacob(2,2)*jacob(3,1)
    tm(6,4) = jacob(3,1)*jacob(1,2) + jacob(3,2)*jacob(1,1)
    tm(4,5) = jacob(1,2)*jacob(2,3) + jacob(1,3)*jacob(2,2)
    tm(5,5) = jacob(2,2)*jacob(3,3) + jacob(2,3)*jacob(3,2)
    tm(6,5) = jacob(3,2)*jacob(1,3) + jacob(3,3)*jacob(1,2)
    tm(4,6) = jacob(1,3)*jacob(2,1) + jacob(1,1)*jacob(2,3)
    tm(5,6) = jacob(2,3)*jacob(3,1) + jacob(2,1)*jacob(3,3)
    tm(6,6) = jacob(3,3)*jacob(1,1) + jacob(3,1)*jacob(1,3)
 
  end subroutine transformation
 
  subroutine get_principal (tensor, eigval, princmatrix)
 
    implicit none
    integer i,j
    real(kind=kreal) :: tensor(1:6)
    real(kind=kreal) :: eigval(3)
    real(kind=kreal) :: princmatrix(3,3)
    real(kind=kreal) :: princnormal(3,3)
    real(kind=kreal) :: tempv(3)
    real(kind=kreal) :: temps
 
    call eigen3(tensor,eigval,princnormal)
 
    if (eigval(1)<eigval(2)) then
      temps=eigval(1)
      eigval(1)=eigval(2)
      eigval(2)=temps
      tempv(:)=princnormal(:,1)
      princnormal(:,1)=princnormal(:,2)
      princnormal(:,2)=tempv(:)
    end if
    if (eigval(1)<eigval(3)) then
      temps=eigval(1)
      eigval(1)=eigval(3)
      eigval(3)=temps
      tempv(:)=princnormal(:,1)
      princnormal(:,1)=princnormal(:,3)
      princnormal(:,3)=tempv(:)
    end if
    if (eigval(2)<eigval(3)) then
      temps=eigval(2)
      eigval(2)=eigval(3)
      eigval(3)=temps
      tempv(:)=princnormal(:,2)
      princnormal(:,2)=princnormal(:,3)
      princnormal(:,3)=tempv(:)
    end if
 
    do j=1,3
      do i=1,3
        princmatrix(i,j) = princnormal(i,j) * eigval(j)
      end do
    end do
 
  end subroutine get_principal
 
  subroutine eigen3d (tensor, eigval, princ)
    implicit none
 
    real(kind=kreal) :: tensor(6)     
    real(kind=kreal) :: eigval(3)     
    real(kind=kreal) :: princ(3,3)   
 
    real(kind=kreal) :: s11, s22, s33, s12, s23, s13, j1, j2, j3
    real(kind=kreal) :: ml,nl
    complex(kind=kreal):: x1,x2,x3
    real(kind=kreal):: rtemp
    integer :: i
    s11 = tensor(1)
    s22 = tensor(2)
    s33 = tensor(3)
    s12 = tensor(4)
    s23 = tensor(5)
    s13 = tensor(6)
 
    ! invariants of stress tensor
    j1 = s11 + s22 + s33
    j2 = -s11*s22 - s22*s33 - s33*s11 + s12**2 + s23**2 + s13**2
    j3 = s11*s22*s33 + 2*s12*s23*s13 - s11*s23**2 - s22*s13**2 - s33*s12**2
    ! Cardano's method
    ! x^3+ ax^2   + bx  +c =0
    ! s^3 - J1*s^2 -J2s -J3 =0
    call cardano(-j1, -j2, -j3, x1, x2, x3)
    eigval(1)= real(x1)
    eigval(2)= real(x2)
    eigval(3)= real(x3)
    if (eigval(1)<eigval(2)) then
      rtemp=eigval(1)
      eigval(1)=eigval(2)
      eigval(2)=rtemp
    end if
    if (eigval(1)<eigval(3)) then
      rtemp=eigval(1)
      eigval(1)=eigval(3)
      eigval(3)=rtemp
    end if
    if (eigval(2)<eigval(3)) then
      rtemp=eigval(2)
      eigval(2)=eigval(3)
      eigval(3)=rtemp
    end if
 
    do i=1,3
      if (eigval(i)/(eigval(1)+eigval(2)+eigval(3)) < 1.0d-10 )then
        eigval(i) = 0.0d0
        princ(i,:) = 0.0d0
        exit
      end if
      ml = ( s23*s13 - s12*(s33-eigval(i)) ) / ( -s23**2 + (s22-eigval(i))*(s33-eigval(i)) )
      nl = ( s12**2 - (s22-eigval(i))*(s11-eigval(i)) ) / ( s12*s23 - s13*(s22-eigval(i)) )
      if (abs(ml) >= huge(ml)) then
        ml=0.0d0
      end if
      if (abs(nl) >= huge(nl)) then
        nl=0.0d0
      end if
      princ(i,1) = eigval(i)/sqrt( 1 + ml**2 + nl**2)
      princ(i,2) = ml * princ(i,1)
      princ(i,3) = nl * princ(i,1)
    end do
 
    write(*,*)
  end subroutine eigen3d
 
  subroutine cardano(a,b,c,x1,x2,x3)
    real(kind=kreal):: a,b,c
    real(kind=kreal):: p,q,d
    complex(kind=kreal):: w
    complex(kind=kreal):: u,v,y
    complex(kind=kreal):: x1,x2,x3
    w = (-1.0d0 + sqrt(dcmplx(-3.0d0)))/2.0d0
    p = -a**2/9.0d0 + b/3.0d0
    q = 2.0d0/2.7d1*a**3 - a*b/3.0d0 + c
    d = q**2 + 4.0d0*p**3
 
    u = ((-dcmplx(q) + sqrt(dcmplx(d)))/2.0d0)**(1.0d0/3.0d0)
 
    if(u.ne.0.0d0) then
      v = -dcmplx(p)/u
      x1 = u + v -dcmplx(a)/3.0d0
      x2 = u*w + v*w**2 -dcmplx(a)/3.0d0
      x3 = u*w**2 + v*w -dcmplx(a)/3.0d0
    else
      y = (-dcmplx(q))**(1.0d0/3.0d0)
      x1 = y -dcmplx(a)/3.0d0
      x2 = y*w -dcmplx(a)/3.0d0
      x3 = y*w**2 -dcmplx(a)/3.0d0
    end if
 
  end subroutine cardano
 
  subroutine rotate_3dvector_by_rodrigues_formula(r,v)
    real(kind=kreal),intent(in)    :: r(3)  
    real(kind=kreal),intent(inout) :: v(3)  
 
    real(kind=kreal) :: rotv(3), rv
    real(kind=kreal) :: cosx, sinc(2) 
    real(kind=kreal) :: x, x2, x4, x6
    real(kind=kreal), parameter :: c0 = 0.5d0
    real(kind=kreal), parameter :: c2 = -4.166666666666666d-002
    real(kind=kreal), parameter :: c4 = 1.388888888888889d-003
    real(kind=kreal), parameter :: c6 = -2.480158730158730d-005
 
    x2 = dot_product(r,r)
    x = dsqrt(x2)
    cosx = dcos(x)
    if( x < 1.d-4 ) then
      x4 = x2*x2
      x6 = x4*x2
      sinc(1) = 1.d0-x2/6.d0+x4/120.d0
      sinc(2) = c0+c2*x2+c4*x4+c6*x6
    else
      sinc(1) = dsin(x)/x
      sinc(2) = (1.d0-cosx)/x2
    endif
 
    ! calc Rot*v
    rv = dot_product(r,v)
    rotv(1:3) = cosx*v(1:3)
    rotv(1:3) = rotv(1:3)+rv*sinc(2)*r(1:3)
    rotv(1) = rotv(1) + (-v(2)*r(3)+v(3)*r(2))*sinc(1)
    rotv(2) = rotv(2) + (-v(3)*r(1)+v(1)*r(3))*sinc(1)
    rotv(3) = rotv(3) + (-v(1)*r(2)+v(2)*r(1))*sinc(1)
    v = rotv
 
  end subroutine
 
  subroutine deriv_general_iso_tensor_func_3d(dpydpx, dydx, eigv, px, py)
    real(kind=kreal), intent(in) :: dpydpx(3,3)
    real(kind=kreal), intent(out) :: dydx(6,6)
    real(kind=kreal), intent(in) :: eigv(3,3)
    real(kind=kreal), intent(in) :: px(3), py(3)
 
    real(kind=kreal) :: x(6)
    real(kind=kreal) :: ep(6,3), dx2dx(6,6)
    integer(kind=kint) :: i, j, m1, m2, m3, ia, ib, ic, ip, jp, k
    real(kind=kreal) :: xmax, dif12, dif23, c1, c2, c3, c4, d1, d2, d3
    real(kind=kreal) :: xa, xc, ya, yc, daa, dac, dca, dcb, dcc
    real(kind=kreal) :: xa_xc, xa_xc2, xa_xc3, ya_yc, xc2
    real(kind=kreal) :: s1, s2, s3, s4, s5, s6
    real(kind=kreal), parameter :: i2(6) = [1.d0, 1.d0, 1.d0, 0.d0, 0.d0, 0.d0]
    real(kind=kreal), parameter :: is(6,6) = reshape([&
        1.d0, 0.d0, 0.d0,  0.d0,   0.d0,   0.d0,&
        0.d0, 1.d0, 0.d0,  0.d0,   0.d0,   0.d0,&
        0.d0, 0.d0, 1.d0,  0.d0,   0.d0,   0.d0,&
        0.d0, 0.d0, 0.d0, 0.5d0,   0.d0,   0.d0,&
        0.d0, 0.d0, 0.d0,  0.d0,  0.5d0,   0.d0,&
        0.d0, 0.d0, 0.d0,  0.d0,   0.d0,  0.5d0],&
        [6, 6])
    real(kind=kreal), parameter :: eps=1.d-6
 
    do i=1,3
      ep(1,i) = eigv(1,i)*eigv(1,i)
      ep(2,i) = eigv(2,i)*eigv(2,i)
      ep(3,i) = eigv(3,i)*eigv(3,i)
      ep(4,i) = eigv(1,i)*eigv(2,i)
      ep(5,i) = eigv(2,i)*eigv(3,i)
      ep(6,i) = eigv(3,i)*eigv(1,i)
    enddo
    x(:) = 0.d0
    do i=1,3
      do j=1,6
        x(j) = x(j)+px(i)*ep(j,i)
      enddo
    enddo
    dx2dx = reshape([2.d0*x(1),      0.d0,      0.d0,             x(4),             0.d0,             x(6),&
        &                 0.d0, 2.d0*x(2),      0.d0,             x(4),             x(5),             0.d0,&
        &                 0.d0,      0.d0, 2.d0*x(3),             0.d0,             x(5),             x(6),&
        &                 x(4),      x(4),      0.d0, (x(1)+x(2))/2.d0,        x(6)/2.d0,        x(5)/2.d0,&
        &                 0.d0,      x(5),      x(6),        x(6)/2.d0, (x(2)+x(3))/2.d0,        x(4)/2.d0,&
        &                 x(6),      0.d0,      x(6),        x(5)/2.d0,        x(4)/2.d0, (x(1)+x(3))/2.d0],&
        [6, 6])
    m1 = maxloc(px, 1)
    m3 = minloc(px, 1)
    if( m1 == m3 ) then
      m1 = 1; m2 = 2; m3 = 3
    else
      m2 = 6 - (m1 + m3)
    endif
    xmax = maxval(abs(px), 1)
    if( xmax < eps ) xmax = 1.d0
    dif12 = abs(px(m1) - px(m2)) / xmax
    dif23 = abs(px(m2) - px(m3)) / xmax
    if( dif12 < eps .and. dif23 < eps ) then
      c1 = dpydpx(1,1)-dpydpx(1,2)
      c2 = dpydpx(1,2)
      do j=1,6
        do i=1,6
          dydx(i,j) = c1*is(i,j)+c2*i2(i)*i2(j)
        enddo
      enddo
    else if( dif12 < eps .or. dif23 < eps ) then
      if( dif12 < eps ) then
        ia = m3; ib = m1; ic = m2
      else
        ia = m1; ib = m2; ic = m3
      endif
      ya = py(ia); yc = py(ic)
      xa = px(ia); xc = px(ic)
      daa = dpydpx(ia,ia)
      dac = dpydpx(ia,ic)
      dca = dpydpx(ic,ia)
      dcb = dpydpx(ic,ib)
      dcc = dpydpx(ic,ic)
      xa_xc = xa-xc
      xa_xc2 = xa_xc*xa_xc
      xa_xc3 = xa_xc2*xa_xc
      ya_yc = ya-yc
      xc2 = xc*xc
      c1 = ya_yc/xa_xc2
      c2 = ya_yc/xa_xc3
      c3 = xc/xa_xc2
      c4 = (xa+xc)/xa_xc
      d1 = dac+dca
      d2 = daa+dcc
      d3 = d1-d2
      s1 = c1+(dcb-dcc)/xa_xc
      s2 = 2.d0*xc*c1+(dcb-dcc)*c4
      s3 = 2.d0*c2+d3/xa_xc2
      s4 = 2.d0*xc*c2+(dac-dcb)/xa_xc+d3*c3
      s5 = 2.d0*xc*c2+(dca-dcb)/xa_xc+d3*c3
      s6 = 2.d0*xc2*c2+(d1*xa-d2*xc)*c3-dcb*c4
      do j=1,6
        do i=1,6
          dydx(i,j) = s1*dx2dx(i,j)-s2*is(i,j)-s3*x(i)*x(j)+s4*x(i)*i2(j)+s5*i2(i)*x(j)-s6*i2(i)*i2(j)
        enddo
      enddo
    else
      do k=1,3
        select case(k)
        case(1)
          ia = 1; ib = 2; ic = 3
        case(2)
          ia = 2; ib = 3; ic = 1
        case(3)
          ia = 3; ib = 1; ic = 2
        end select
        c1 = py(ia)/((px(ia)-px(ib))*(px(ia)-px(ic)))
        c2 = px(ib)+px(ic)
        c3 = (px(ia)-px(ib))+(px(ia)-px(ic))
        c4 = px(ib)-px(ic)
        do j=1,6
          do i=1,6
            dydx(i,j) = dydx(i,j)+c1*(dx2dx(i,j)-c2*is(i,j)-c3*ep(i,ia)*ep(j,ia)-c4*(ep(i,ib)*ep(j,ib)-ep(i,ic)*ep(j,ic)))
          enddo
        enddo
      enddo
      do jp=1,3
        do ip=1,3
          c1 = dpydpx(ip,jp)
          do j=1,6
            do i=1,6
              dydx(i,j) = dydx(i,j)+c1*ep(i,ip)*ep(j,jp)
            enddo
          enddo
        enddo
      enddo
    endif
  end subroutine deriv_general_iso_tensor_func_3d
 
end module