kinematic.f90.svn-base

! ------------------------------------------------------------------------
!>  \brief Module for Kinematic Matrix
! ------------------------------------------------------------------------
MODULE kinematic
USE system_constants    !> Global constants (for size of matrix)
USE dense_matrix_def    !> Dense Matrix data type
USE field_data          !> Field Data module

IMPLICIT NONE
PRIVATE



! ************************************************************************
! EXPORTS
! ************************************************************************


!> Exported interfaces
PUBLIC :: bmatrix



! ************************************************************************
! LOCAL CONSTANTS
! ************************************************************************

! none


! ************************************************************************
! DATA TYPES
! ************************************************************************

! none



! ************************************************************************
! INTERFACES
! ************************************************************************


!> \brief Interface to kinematic matrix access program
INTERFACE bmatrix
  MODULE PROCEDURE bmatrix_3_noded_small_strain_
END INTERFACE bmatrix



CONTAINS



! ************************************************************************
! ACCESS PROGRAMS
! ************************************************************************


! ------------------------------------------------------------------------
!>  \brief Kinematic matrix (3-noded, small strain)
!!
!!  \param  iel   Element number
!!  \param  Bmat  Kinematic matrix
!!
!!  This routine builds the kinematic matrix for a 3-noded triangular
!!  element considering small strain. This is sometimes referred to in the
!!  finite element literature as the B matrix, which is defined as:
!!
!!        Bmat = L*N
!!
!!  where L is the associated linear differential operator and N is the
!!  matrix of shape functions. Using the mathematical definition of
!!  strain:
!!
!!        L = [ d/dx     0   ]
!!            [    0   d/dy  ]
!!            [ d/dy   d/dx  ]
!!
!!  where d/dx and d/dy are partial derivatives with respect to the
!!  coordinate axes. Combining this with the shape functions based on
!!  area coordinates for a 3-noded triangular element, the kinematic
!!  matrix is:
!!
!!        Bmat = ( 1/(2*A) ) *  [ b1    0   b2    0   b3    0  ]
!!                              [  0   c1    0   c2    0   c3  ]
!!                              [ c1   b1   c2   b2   c3   b3  ]
!!
!!  where bi and ci are defined as:
!!
!!        bi = yj - yk      for {i,j,k} \in {{1,2,3},{2,3,1},{3,1,2}}
!!        ci = xk - xj
!!
!!  and (xi,yi) are the nodal coordinates of the element.
! ------------------------------------------------------------------------
SUBROUTINE bmatrix_3_noded_small_strain_ (iel, Bmat)
  INTEGER, INTENT(IN) :: iel
  TYPE(matrixT), INTENT(INOUT) :: Bmat
  DOUBLE PRECISION, DIMENSION(NNODEL) :: x,y  !> for nodal coords
  DOUBLE PRECISION, DIMENSION(NNODEL) :: b,c  !> for coefficients
  DOUBLE PRECISION :: A       !> area of element
  DOUBLE PRECISION :: coef    !> area coefficient
  INTEGER :: i    !> loop variable

  !> get area of element
  A = fld_volElem(iel)
  coef = 1.d0/(2.d0*A)

  !> get nodal coords
  DO i = 1,NNODEL
    x(i) = fld_getCoord( fld_getConnect(iel,i), 1 )
    y(i) = fld_getCoord( fld_getConnect(iel,i), 2 )
  END DO

  !> compute coefficients
  b(1) = y(2)-y(3)
  b(2) = y(3)-y(1)
  b(3) = y(1)-y(2)
  b = coef*b

  c(1) = x(3)-x(2)
  c(2) = x(1)-x(3)
  c(3) = x(2)-x(1)
  c = coef*c

  !> initialize matrix
  CALL dm_init(Bmat, NTNS,NDIM*NNODEL)

  !> set matrix contents
  Bmat%dat(1,1:NDIM*NNODEL:2) = b(1:NNODEL)
  Bmat%dat(2,2:NDIM*NNODEL:2) = c(1:NNODEL)
  Bmat%dat(3,1:NDIM*NNODEL:2) = c(1:NNODEL)
  Bmat%dat(3,2:NDIM*NNODEL:2) = b(1:NNODEL)

END SUBROUTINE bmatrix_3_noded_small_strain_



! ************************************************************************
! LOCAL FUNCTIONS
! ************************************************************************

! none

END MODULE kinematic