StiffnessMatrixClass.f90 Source File

module StiffnessMatrixClass
   use ShapeFunctionClass
   use MeshClass
   use StressClass
   use StrainClass

   implicit none

   ! Provides a set of Element-by-element stiffness matrix & RHS vector of internal force
   ! for all types of elements, all Strain Theory, and all constitutive models.

   type :: StiffnessMatrix_
      real(8), allocatable     :: Amat(:, :), bvec(:)
      type(ShapeFunction_)    :: ShapeFunc
      type(Stress_), pointer   :: Stress
      type(Strain_), pointer   :: Strain

      integer :: TheoryID, ElemID, GpID

      character*40 :: StrainTheory

      ! Please input one of following keyphrases

      ! Finite_Elasticity
      ! Finite_ElastoPlasticity
      ! Infinitesimal_Elasticity
      ! Infinitesimal_ElastoPlasticity
      ! Small_strain
   contains
      procedure, public :: init => initStiffnessMatrix
      procedure, public :: update => updateStiffnessMatrix
   end type

contains

! ######################################
   subroutine initStiffnessMatrix(obj, FEMDomain, Stress, Strain, withInit)
      class(StiffnessMatrix_), intent(inout)::obj
      !character(*),intent(in) :: StrainTheory
      class(FEMDomain_), intent(in) :: FEMDomain
      class(Stress_), target, intent(in) :: Stress
      class(Strain_), target, intent(in) :: Strain
      logical, optional, intent(in) :: withInit
      integer :: NumOfDim, NumOfNodePerElem

      obj%Stress => Stress
      obj%Strain => Strain
      ! create Stress/Strain measure objects
      !obj%StrainTheory = StrainTheory
      !call obj%Stress%init(StrainTheory)
      !call obj%Stress%init(StrainTheory)

      ! create a ShapeFunction object
      NumOfDim = size(FEMDomain%Mesh%NodCoord, 2)
      NumOfNodePerElem = size(FEMDomain%Mesh%ElemNod, 2)
      call obj%ShapeFunc%getType(NumOfDim=NumOfDim, NumOfNodePerElem=NumOfNodePerElem)
      call obj%ShapeFunc%setType()

      obj%ElemID = 1
      obj%GpID = 1

   end subroutine
! ######################################

! ######################################
   subroutine updateStiffnessMatrix(obj, Mesh, ElemID, MeshID)
      class(StiffnessMatrixClass_), intent(inout) :: obj
      class(Mesh_), intent(in) :: Mesh
      integer, optional, intent(in) :: ElemID, MeshID
      integer :: ElemID_out, MeshID_out

      ! get element ID and Gauss-point-id
      if (present(ElemID)) then
         ElemID_out = ElemID
      else
         ElemID_out = obj%ElemID
      end if
      if (present(GpID)) then
         GpID_out = GpID
      else
         GpID_out = obj%GpID
      end if

      ! get all objects related to a shape function
      call obj%ShapeFunc%update(ElemType=obj%ShapeFunc%ElemType, &
                                NodCoord=Mesh%NodCoordInit, NodCoord=Mesh%NodCoord_n, &
                                ElemNod=Mesh%ElemNod, ElemID=ElemID_out, GpID=GpID_out)

      ! create strain tensor

      ! create stress tensor

      ! create stiffness matrix

      ! update element ID and Gauss-point-id
      obj%GpID = obj%GpID + 1
      if (obj%GpID > size(obj%ShapeFunc%NumOfGp)) then
         ! go to next element
         obj%GpID = 1
         obj%ElemID = obj%ElemID + 1
      else
         ! go to next Gauss point in the same element
         obj%GpID = obj%GpID + 1
      end if
   end subroutine
! ######################################

end module