v3d10_thermalExp.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE v3d10_thermalexp(coor,matcstet,lmcstet,R_in,ci, &
     s11l,s22l,s33l,&
     numnp,nstart,nend,numcstet,numat_vol, coeffexp,temperature,temperature0)

!!****f* Rocfrac/Rocfrac/Source/v3d10_thermalExp.f90
!!
!!  NAME
!!    v3d10_thermalExp
!!
!!  FUNCTION
!!
!!    Computes the internal force for a small deformation, 
!!    linear elastic 10-node tetrahedral with thermal expansion
!!
!!  INPUTS
!!
!!   NumNP    -- Number of nodes
!!   numcstet -- Number of elements
!!   Coor     -- number of coordinates
!!   Matcstet  -- Material id
!!   d -- nodal displacement
!!   ci -- compliance matrix
!! 
!!       stored in the follow format:
!!
!!            [ ci(1) ci(2) ci(4)                   ]
!!            [       ci(3) ci(5)                   ]
!!            [             ci(6)                   ]
!!            [                   ci(7)             ]
!!            [                         ci(8)       ]
!!            [                               ci(9) ]
!!
!!   CoeffExp -- coeffecient of thermal expansion
!!   Temperature - nodal temperature
!!   Temperature0 - reference temperature
!!    
!!
!!   lmcstet  -- Nodal connectivity
!!   nstart, nend -- element beginning and end loop counter
!!   numat_vol -- number of materials
!!
!!  OUTPUT
!!
!!    Rin -- internal force vector
!!    S11l,S22l,S33l,S12l,S23l,S13l -- Stresses at gauss points
!!
!!****

!----- 4 Guass points 
  IMPLICIT NONE
!-----Global variables
  INTEGER :: numnp          ! number of nodes
  INTEGER :: numat_vol      ! number of volumetric materials
  INTEGER :: numcstet       ! number of LSTets
!--   coordinate array
  REAL*8, DIMENSION(1:3,1:numnp) :: coor
!--   elastic stiffness consts
  REAL*8, DIMENSION(1:9,1:numat_vol) :: ci

  REAL*8, DIMENSION(1:numat_vol) :: coeffexp

!--   internal force
  REAL*8, DIMENSION(1:3*numnp) :: r_in
!--   CSTet stress
  REAL*8, DIMENSION(1:4,1:numcstet) :: s11l, s22l, s33l

  REAL*8, DIMENSION(1:4) :: strssth11, strssth22, strssth33


!--   connectivity table for CSTet  
  INTEGER, DIMENSION(1:10,1:numcstet) :: lmcstet
!--   mat number for CSTet element
  INTEGER, DIMENSION(1:numcstet) :: matcstet
!---- Local variables
!--   node numbers
  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
!--   6*volume and inverse of 6*volume      
  REAL*8 :: vx6, vx6inv
!--   spacial derivatives
  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
  REAL*8 :: b11,b12,b13,b14,b15,b16,b17,b18,b19,b20
  REAL*8 :: b21,b22,b23,b24,b25,b26,b27,b28,b29,b30
!--   strains
  REAL*8 :: e11,e22,e33,e12,e23,e13
!--   coordinate holding variable
  REAL*8 :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
  REAL*8 :: y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
  REAL*8 :: z1,z2,z3,z4,z5,z6,z7,z8,z9,z10
!--   dummy and counters
  INTEGER :: i,j,nstart,nend
  REAL*8 :: aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12
!--   partial internal force
  REAL*8 :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,r17,r18
  REAL*8 :: r19,r20,r21,r22,r23,r24,r25,r26,r27,r28,r29,r30
  REAL*8 :: g1, g2, g3, g4
  REAL*8 :: xn1, xn2, xn3, xn4
!--  Coordinate subtractions
  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
!--  
  REAL*8 :: c11, c21, c31
  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8,k1n9,k1n10
  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8,k2n9,k2n10
  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8,k3n9,k3n10

  REAL*8, DIMENSION(1:NumNP) :: temperature
  REAL*8 :: temperaturegauss,temperature0


  REAL*8,DIMENSION(1:4,1:4) :: gaussintpt = reshape( &
       (/0.58541020d0,0.13819660d0,0.13819660d0,0.13819660d0, &
       0.13819660d0,0.58541020d0,0.13819660d0,0.13819660d0, &
       0.13819660d0,0.13819660d0,0.58541020d0,0.13819660d0, &
       0.13819660d0,0.13819660d0,0.13819660d0,0.58541020d0/),(/4,4/) )
  
  INTEGER :: igpt

  DO i = nstart, nend

     j   = matcstet(i)
     
     n1  = lmcstet(1,i)
     n2  = lmcstet(2,i)
     n3  = lmcstet(3,i)
     n4  = lmcstet(4,i)
     n5  = lmcstet(5,i)
     n6  = lmcstet(6,i)
     n7  = lmcstet(7,i)
     n8  = lmcstet(8,i)
     n9  = lmcstet(9,i)
     n10 = lmcstet(10,i)
     
     k3n1  = 3*n1
     k3n2  = 3*n2
     k3n3  = 3*n3
     k3n4  = 3*n4
     k3n5  = 3*n5
     k3n6  = 3*n6
     k3n7  = 3*n7
     k3n8  = 3*n8
     k3n9  = 3*n9
     k3n10 = 3*n10
     
     k2n1  = k3n1  - 1
     k2n2  = k3n2  - 1
     k2n3  = k3n3  - 1
     k2n4  = k3n4  - 1
     k2n5  = k3n5  - 1
     k2n6  = k3n6  - 1
     k2n7  = k3n7  - 1
     k2n8  = k3n8  - 1
     k2n9  = k3n9  - 1
     k2n10 = k3n10 - 1
     
     k1n1  = k3n1  - 2
     k1n2  = k3n2  - 2
     k1n3  = k3n3  - 2
     k1n4  = k3n4  - 2 
     k1n5  = k3n5  - 2
     k1n6  = k3n6  - 2
     k1n7  = k3n7  - 2
     k1n8  = k3n8  - 2 
     k1n9  = k3n9  - 2
     k1n10 = k3n10 - 2
                                ! k#n# dummy variables replaces:
     x1 = coor(1,n1)
     x2 = coor(1,n2)
     x3 = coor(1,n3)
     x4 = coor(1,n4)
     y1 = coor(2,n1)
     y2 = coor(2,n2)
     y3 = coor(2,n3)
     y4 = coor(2,n4)
     z1 = coor(3,n1)
     z2 = coor(3,n2)
     z3 = coor(3,n3)
     z4 = coor(3,n4)
     
     x14 = x1 - x4
     x24 = x2 - x4
     x34 = x3 - x4
     y14 = y1 - y4
     y24 = y2 - y4
     y34 = y3 - y4
     z14 = z1 - z4
     z24 = z2 - z4
     z34 = z3 - z4
     
     c11 =    y24*z34 - z24*y34
     c21 = -( x24*z34 - z24*x34 )
     c31 =    x24*y34 - y24*x34

     vx6 = -( x14*c11 + y14*c21 + z14*c31 )

!        Vx6 = x2*y3*z4 - x2*y4*z3 - y2*x3*z4 + y2*x4*z3 + z2*x3*y4 
!     $        - z2*x4*y3 - x1*y3*z4 + x1*y4*z3 + x1*y2*z4 - x1*y2*z3
!     $        - x1*z2*y4 + x1*z2*y3 + y1*x3*z4 - y1*x4*z3 - y1*x2*z4 
!     $        + y1*x2*z3 + y1*z2*x4 - y1*z2*x3 - z1*x3*y4 + z1*x4*y3
!     $        + z1*x2*y4 - z1*x2*y3 - z1*y2*x4 + z1*y2*x3

     vx6inv = 1.d0 / vx6
     
     aux1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3)
     aux2 =  (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3)
     aux3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3)
     aux4 =  (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3)
     aux5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3)
     aux6 =  (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3)
     aux7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2)
     aux8 =  (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2)
     aux9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2)
     aux10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2)
     aux11 =-(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2)
     aux12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2)

     r1 = 0.
     r2 = 0.
     r3 = 0.
     r4 = 0.
     r5 = 0.
     r6 = 0.
     r7 = 0.
     r8 = 0.
     r9 = 0.
     r10 = 0.
     r11 = 0.
     r12 = 0.
     r13 = 0.
     r14 = 0.
     r15 = 0.
     r16 = 0.
     r17 = 0.
     r18 = 0.
     r19 = 0.
     r20 = 0.
     r21 = 0.
     r22 = 0.
     r23 = 0.
     r24 = 0.
     r25 = 0.
     r26 = 0.
     r27 = 0.
     r28 = 0.
     r29 = 0.
     r30 = 0.

     DO igpt = 1, 4

        g1 = gaussintpt(igpt,1)
        g2 = gaussintpt(igpt,2)
        g3 = gaussintpt(igpt,3)
        g4 = gaussintpt(igpt,4)

        xn1 = (4.d0*g1-1.d0)   ! derivative of shape function
        xn2 = (4.d0*g2-1.d0)   ! dN_i/dzeta_i
        xn3 = (4.d0*g3-1.d0)
        xn4 = (4.d0*g4-1.d0)
!        xN5 = 4.d0*g1*g2
!        xN6 = 4.d0*g2*g3
!        xN7 = 4.d0*g3*g1
!        xN8 = 4.d0*g1*g4
!        xN9 = 4.d0*g2*g4
!        xN10= 4.d0*g3*g4

        b1 = aux1*xn1
        b2 = aux2*xn1
        b3 = aux3*xn1
        b4 = aux4*xn2
        b5 = aux5*xn2
        b6 = aux6*xn2
        b7 = aux7*xn3
        b8 = aux8*xn3
        b9 = aux9*xn3
        b10 =  aux10*xn4
        b11 =  aux11*xn4
        b12 =  aux12*xn4
        
        b13 =  4.d0*(g2*aux1 + g1*aux4)
        b14 =  4.d0*(g2*aux2 + g1*aux5)
        b15 =  4.d0*(g2*aux3 + g1*aux6)
        
        b16 =  4.d0*(g3*aux4 + g2*aux7)
        b17 =  4.d0*(g3*aux5 + g2*aux8)
        b18 =  4.d0*(g3*aux6 + g2*aux9)
        
        b19 =  4.d0*(g1*aux7 + g3*aux1)
        b20 =  4.d0*(g1*aux8 + g3*aux2)
        b21 =  4.d0*(g1*aux9 + g3*aux3)
        
        b22 =  4.d0*(g4*aux1 + g1*aux10)
        b23 =  4.d0*(g4*aux2 + g1*aux11)
        b24 =  4.d0*(g4*aux3 + g1*aux12)
        
        b25 =  4.d0*(g4*aux4 + g2*aux10)
        b26 =  4.d0*(g4*aux5 + g2*aux11)
        b27 =  4.d0*(g4*aux6 + g2*aux12)
        
        b28 =  4.d0*(g4*aux7 + g3*aux10)
        b29 =  4.d0*(g4*aux8 + g3*aux11)
        b30 =  4.d0*(g4*aux9 + g3*aux12)
        

! interpolate Temperature from nodes to Gauss Point

        temperaturegauss = g1*(2.d0*g1-1.d0)*temperature(n1) +                 &
             g2*(2.d0*g2-1.d0)*temperature(n2) +                               &
             g3*(2.d0*g3-1.d0)*temperature(n3) +                               &
             g4*(2.d0*g4-1.d0)*temperature(n4) +                               &
             4.d0*g1*g2*temperature(n5) + 4.d0*g2*g3*temperature(n6) +         &
             4.d0*g3*g1*temperature(n7) + 4.d0*g1*g4*temperature(n8) +         &
             4.d0*g2*g4*temperature(n9) + 4.d0*g3*g4*temperature(n10)

        strssth11(igpt) = -ci(1,j)*coeffexp(j)*( temperaturegauss - temperature0 ) !* Vx6inv
        strssth22(igpt) = -ci(3,j)*coeffexp(j)*( temperaturegauss - temperature0 ) !* Vx6inv
        strssth33(igpt) = -ci(6,j)*coeffexp(j)*( temperaturegauss - temperature0 ) !* Vx6inv


        r1 = r1 + strssth11(igpt)*b1
        r2 = r2 + strssth22(igpt)*b2
        r3 = r3 + strssth33(igpt)*b3
     
        r4 = r4 + strssth11(igpt)*b4
        r5 = r5 + strssth22(igpt)*b5
        r6 = r6 + strssth33(igpt)*b6
        
        r7 = r7 + strssth11(igpt)*b7 
        r8 = r8 + strssth22(igpt)*b8 
        r9 = r9 + strssth33(igpt)*b9 
        
        r10 = r10 + strssth11(igpt)*b10 
        r11 = r11 + strssth22(igpt)*b11 
        r12 = r12 + strssth33(igpt)*b12 
        
        r13 = r13 + strssth11(igpt)*b13 
        r14 = r14 + strssth22(igpt)*b14 
        r15 = r15 + strssth33(igpt)*b15 
        
        r16 = r16 + strssth11(igpt)*b16 
        r17 = r17 + strssth22(igpt)*b17 
        r18 = r18 + strssth33(igpt)*b18 
        
        r19 = r19 + strssth11(igpt)*b19 
        r20 = r20 + strssth22(igpt)*b20 
        r21 = r21 + strssth33(igpt)*b21 
        
        r22 = r22 + strssth11(igpt)*b22 
        r23 = r23 + strssth22(igpt)*b23 
        r24 = r24 + strssth33(igpt)*b24 
        
        r25 = r25 + strssth11(igpt)*b25 
        r26 = r26 + strssth22(igpt)*b26 
        r27 = r27 + strssth33(igpt)*b27 
        
        r28 = r28 + strssth11(igpt)*b28 
        r29 = r29 + strssth22(igpt)*b29 
        r30 = r30 + strssth33(igpt)*b30 

        s11l(igpt,i) = s11l(igpt,i) + strssth11(igpt)
        s22l(igpt,i) = s22l(igpt,i) + strssth22(igpt)
        s33l(igpt,i) = s33l(igpt,i) + strssth33(igpt)

     ENDDO

!     PRINT*,i, S11l(1:4,i)

! Wi (i.e. weight) for 4 guass integration is 1/4         

! Wi * 1/6 because the volume of a reference tetrahedra in 
! volume coordinates is 1/6

! ASSEMBLE THE INTERNAL FORCE VECTOR
!
! local node 1
     r_in(k1n1)  = r_in(k1n1)  - r1*0.04166666666666667d0
     r_in(k2n1)  = r_in(k2n1)  - r2*0.04166666666666667d0
     r_in(k3n1)  = r_in(k3n1)  - r3*0.04166666666666667d0
! local node 2
     r_in(k1n2)  = r_in(k1n2)  - r4*0.04166666666666667d0
     r_in(k2n2)  = r_in(k2n2)  - r5*0.04166666666666667d0
     r_in(k3n2)  = r_in(k3n2)  - r6*0.04166666666666667d0
! local node 3
     r_in(k1n3)  = r_in(k1n3)  - r7*0.04166666666666667d0
     r_in(k2n3)  = r_in(k2n3)  - r8*0.04166666666666667d0
     r_in(k3n3)  = r_in(k3n3)  - r9*0.04166666666666667d0
! local node 4
     r_in(k1n4)  = r_in(k1n4)  - r10*0.04166666666666667d0
     r_in(k2n4)  = r_in(k2n4)  - r11*0.04166666666666667d0
     r_in(k3n4)  = r_in(k3n4)  - r12*0.04166666666666667d0
! local node 5
     r_in(k1n5)  = r_in(k1n5)  - r13*0.04166666666666667d0
     r_in(k2n5)  = r_in(k2n5)  - r14*0.04166666666666667d0
     r_in(k3n5)  = r_in(k3n5)  - r15*0.04166666666666667d0
! local node 6
     r_in(k1n6)  = r_in(k1n6)  - r16*0.04166666666666667d0
     r_in(k2n6)  = r_in(k2n6)  - r17*0.04166666666666667d0
     r_in(k3n6)  = r_in(k3n6)  - r18*0.04166666666666667d0
! local node 7
     r_in(k1n7)  = r_in(k1n7)  - r19*0.04166666666666667d0
     r_in(k2n7)  = r_in(k2n7)  - r20*0.04166666666666667d0
     r_in(k3n7)  = r_in(k3n7)  - r21*0.04166666666666667d0
! local node 8
     r_in(k1n8)  = r_in(k1n8)  - r22*0.04166666666666667d0
     r_in(k2n8)  = r_in(k2n8)  - r23*0.04166666666666667d0
     r_in(k3n8)  = r_in(k3n8)  - r24*0.04166666666666667d0
! local node 9
     r_in(k1n9)  = r_in(k1n9)  - r25*0.04166666666666667d0
     r_in(k2n9)  = r_in(k2n9)  - r26*0.04166666666666667d0
     r_in(k3n9)  = r_in(k3n9)  - r27*0.04166666666666667d0
! local node 10
     r_in(k1n10) = r_in(k1n10) - r28*0.04166666666666667d0
     r_in(k2n10) = r_in(k2n10) - r29*0.04166666666666667d0
     r_in(k3n10) = r_in(k3n10) - r30*0.04166666666666667d0

  ENDDO
  
  RETURN
END SUBROUTINE v3d10_thermalexp