jacobview_v.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    jacobiew_v       ../engine/source/materials/mat_share/jacobview_v.F
!||--- called by ------------------------------------------------------
!||    fail_inievo_ib   ../engine/source/materials/fail/inievo/fail_inievo_ib.F90
!||    fail_inievo_s    ../engine/source/materials/fail/inievo/fail_inievo_s.F
!||====================================================================
        SUBROUTINE  jacobiew_v(DIM1,DIM2,A,EW,EV,NROT)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
        INTEGER, INTENT(in) :: DIM1
        INTEGER, INTENT(in) :: DIM2
        INTEGER, DIMENSION(DIM1) :: NROT
        my_real, DIMENSION(DIM1,DIM2), intent(inout) :: ew 
        my_real, DIMENSION(DIM1,DIM2,DIM2), intent(inout) :: a,ev
 
        my_real, DIMENSION(:,:), ALLOCATABLE :: b,z
        INTEGER IZ,IS,ITER,J,IJK,I
        INTEGER :: NINDX
        INTEGER, DIMENSION(DIM1) :: INDX
        my_real, DIMENSION(DIM1) :: sumrs,eps
        my_real :: g,h,t,c,s,tau,theta      
C----------------------------------------------------------------
 
        ! ---------------------
        ALLOCATE( b(dim1,dim2) )
        ALLOCATE( z(dim1,dim2) )
        DO iz=1,dim2
            DO is=1,dim2
                IF(iz>is) THEN 
                    DO i=1,dim1
                        a(i,is,iz) = a(i,iz,is) 
                    ENDDO
                ENDIF
                ev(1:dim1,iz,is)=zero
            ENDDO
            b(1:dim1,iz)=a(1:dim1,iz,iz)
            ew(1:dim1,iz)=b(1:dim1,iz)
            z(1:dim1,iz)=0.
            ev(1:dim1,iz,iz)=one
        ENDDO
        nrot(1:dim1)=0
        ! ---------------------
 
 
        ! ---------------------
        DO iter = 1,50
            sumrs(1:dim1) = zero
            ! ---------------------  
            ! diagonal sum          
            DO iz=1,dim2-1
                DO is=iz+1,dim2
                    sumrs(1:dim1)=sumrs(1:dim1)+abs(a(1:dim1,iz,is))
                ENDDO
            ENDDO
            ! ---------------------
            ! check if the computation is mandatory : sumrs(i) = 0 --> skip the element
            nindx = 0
            DO i=1,dim1
                IF (sumrs(i)/=zero) THEN
                    nindx = nindx + 1
                    indx(nindx) = i
                ENDIF
            ENDDO
            ! ---------------------
            IF (iter > 4)   THEN
                eps(1:dim1) = zero
            ELSE
                eps(1:dim1) = one_fifth*sumrs(1:dim1)/dim2**2
            ENDIF
            ! ---------------------
            DO iz=1,dim2-1
                DO is=iz+1,dim2
#include "vectorize.inc"                       
                    DO ijk=1,nindx
                        i = indx(ijk)
                        g = 100. * abs(a(i,iz,is))
                        IF ( iter>4 .AND. abs(ew(i,iz))+g==abs(ew(i,iz))
     &                          .AND. abs(ew(i,is))+g==abs(ew(i,is))) THEN
                            a(i,iz,is)=zero
                        ELSEIF (abs(a(i,iz,is)) > eps(i)) THEN
                            h = ew(i,is)-ew(i,iz)
                            IF (abs(h)+g==abs(h)) THEN
                                t = a(i,iz,is)/h
                            ELSE
                                theta = half*h/a(i,iz,is)
                                t=one/(abs(theta)+sqrt(one+theta**2))
                                IF (theta < zero) t=-t
                            ENDIF
                            c=one/sqrt(one+t**2)
                            s=t*c
                            tau=s/(one+c)
                            h=t*a(i,iz,is)
                            z(i,iz)=z(i,iz)-h
                            z(i,is)=z(i,is)+h
                            ew(i,iz)=ew(i,iz)-h
                            ew(i,is)=ew(i,is)+h
                            a(i,iz,is)=zero
                            DO j=1,iz-1
                                g=a(i,j,iz)
                                h=a(i,j,is)
                                a(i,j,iz)=g-s*(h+g*tau)
                                a(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            DO j=iz+1,is-1
                                g=a(i,iz,j)
                                h=a(i,j,is)
                                a(i,iz,j)=g-s*(h+g*tau)
                                a(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            DO j=is+1,dim2
                                g=a(i,iz,j)
                                h=a(i,is,j)
                                a(i,iz,j)=g-s*(h+g*tau)
                                a(i,is,j)=h+s*(g-h*tau)
                            ENDDO
                            DO j=1,dim2
                                g=ev(i,j,iz)
                                h=ev(i,j,is)
                                ev(i,j,iz)=g-s*(h+g*tau)
                                ev(i,j,is)=h+s*(g-h*tau)
                            ENDDO
                            nrot(i)=nrot(i)+1                        
                        ENDIF
                    ENDDO
                ENDDO
            ENDDO
            ! ---------------------
            ! update b/ew/z
            DO iz=1,dim2
#include "vectorize.inc"                    
                DO ijk=1,nindx
                    i = indx(ijk)
                    b(i,iz)=b(i,iz)+z(i,iz)
                    ew(i,iz)=b(i,iz)
                    z(i,iz)=zero
                ENDDO
            ENDDO
            ! ---------------------
        ENDDO
 
        DEALLOCATE( b )
        DEALLOCATE( z )
 
        RETURN
        END SUBROUTINE jacobiew_v