module m_MEF90_Materials_Types #include "petsc/finclude/petsc.h" use m_MEF90_Parameters use m_MEF90_Utils use m_MEF90_LinAlg implicit none(type, external) type MEF90HookesLaw2D type(Tens4OS2D) :: fullTensor type(Tens4OS3D) :: fullTensorLocal, fullTensor3D PetscReal :: lambda, mu, YoungsModulus, PoissonRatio, BulkModulus PetscEnum :: type PetscBool :: isPlaneStress end type MEF90HookesLaw2D type MEF90HookesLaw3D type(Tens4OS3D) :: fullTensor, fullTensorLocal PetscReal :: lambda, mu, YoungsModulus, PoissonRatio, BulkModulus PetscEnum :: type #if (PETSC_SIZEOF_INT == 4) ! With 4-byte integers, this declared type is 4-bytes shy of being ! 8-byte aligned, which can be problematic for arrays of this type. ! The MEF90HookesLaw2D has an extra PetscBool that keeps it 8-byte ! aligned, so we'll mimic that: PetscBool :: padding = PETSC_FALSE #endif end type MEF90HookesLaw3D type MEF90RotationMatrix3D ! The Bunge (passive) convention is used. The rotation matrix transforms a vector from the global frame to the local frame. ! - rotation of a vector: [V_local]_i = [R]_ij . [V_global]_j ! - rotation of a 2nd order tensor: [M_local]_ij = [R]_ik . [M_global]_kl . [R^T]_lj ! - rotation of a fourth order tensor: [C_local]_ijkl = [R]_ip . [R]_jq . [R]_kr . [R]_ls . [C_global]_pqrs type(MAT3D) :: fullTensor PetscReal :: phi1, Phi, phi2 type(Vect3D) :: V1, V2, V3 PetscBool :: fromEuler end type MEF90RotationMatrix3D enum, bind(c) enumerator :: MEF90HookesLawTypeFull = 0, & MEF90HookesLawTypeIsotropic end enum end module m_MEF90_Materials_Types module m_MEF90_Materials #include "petsc/finclude/petsc.h" use m_MEF90_Parameters use m_MEF90_Utils use m_MEF90_LinAlg use m_MEF90_Elements use m_MEF90_Ctx use m_MEF90_DMPlex use m_MEF90_Materials_Types implicit none(type) interface operator(+) module procedure MEF90HookesLaw2DSum, MEF90HookesLaw3DSum end interface interface operator(-) module procedure MEF90HookesLaw2DDiff, MEF90HookesLaw3DDiff end interface interface operator(*) module procedure MEF90HookesLaw2DXMatS2D, MEF90HookesLaw3DXMatS3D, MEF90HookesLaw2DXMat2D, MEF90HookesLaw3DXMat3D, ScalarXMEF90HookesLaw2D, ScalarXMEF90HookesLaw3D end interface PetscSizeT, protected :: sizeofMEF90HookesLaw2D PetscSizeT, protected :: sizeofMEF90HookesLaw3D character(len=MEF90MXSTRLEN), dimension(5), protected :: MEF90HookesLawTypeList contains #undef __FUNCT__ #define __FUNCT__ "MEF90MaterialsInitialize_Private" !!! author: Blaise Bourdin (2013-2014, bourdin@lsu.edu) !!! !!! MEF90MaterialsInitialize_Private: !!! subroutine MEF90MaterialsInitialize_Private(ierr) PetscErrorCode, intent(OUT) :: ierr type(MEF90HookesLaw2D), target :: HookesLaw2D type(MEF90HookesLaw3D), target :: HookesLaw3D character(len=1), pointer :: dummychar(:) PetscSizeT :: sizeofchar PetscCall(PetscDataTypeGetSize(PETSC_CHAR, sizeofchar, ierr)) sizeofMEF90HookesLaw2D = size(transfer(HookesLaw2D, dummychar)) * sizeofchar sizeofMEF90HookesLaw3D = size(transfer(HookesLaw3D, dummychar)) * sizeofchar MEF90HookesLawTypeList(1) = 'Full' MEF90HookesLawTypeList(2) = 'Isotropic' MEF90HookesLawTypeList(3) = 'MEF90HookesLawTypeList' MEF90HookesLawTypeList(4) = '_MEF90HookesLawTypeList' MEF90HookesLawTypeList(5) = '' end subroutine MEF90MaterialsInitialize_Private !!! Subroutine generating various types of Hooke's laws #undef __FUNCT__ #define __FUNCT__ "MEF90HookeLawIsoLambdaMu2D" subroutine MEF90HookeLawIsoLambdaMu2D(A, lambda, mu) type(Tens4OS2D), intent(OUT) :: A PetscReal, intent(IN) :: lambda, mu A = 0.0_kr A%XXXX = lambda + 2.0_kr * mu A%XXYY = lambda A%XYXY = mu A%YYYY = lambda + 2.0_kr * mu end subroutine MEF90HookeLawIsoLambdaMu2D #undef __FUNCT__ #define __FUNCT__ "MEF90HookeLawIsoEnu2DPlaneStress" subroutine MEF90HookeLawIsoEnu2DPlaneStress(A, E, nu) PetscReal, intent(IN) :: E, nu type(Tens4OS2D), intent(OUT) :: A PetscReal :: Lambda, mu lambda = E * nu / (1.0_kr - nu**2) mu = E / (1.0_kr + nu)*.5_kr A = 0.0_kr A%XXXX = lambda + 2.0_kr * mu A%XXYY = lambda A%XYXY = mu A%YYYY = lambda + 2.0_kr * mu end subroutine MEF90HookeLawIsoEnu2DPlaneStress #undef __FUNCT__ #define __FUNCT__ "MEF90HookeLawIsoEnu2DPlaneStrain" subroutine MEF90HookeLawIsoEnu2DPlaneStrain(A, E, nu) PetscReal, intent(IN) :: E, nu type(Tens4OS2D), intent(OUT) :: A PetscReal :: Lambda, mu lambda = E * nu / (1.0_kr + nu) / (1.0_kr - 2.0_kr * nu) mu = E / (1.0_kr + nu)*.5_kr A = 0.0_kr A%XXXX = lambda + 2.0_kr * mu A%XXYY = lambda A%XYXY = mu A%YYYY = lambda + 2.0_kr * mu end subroutine MEF90HookeLawIsoEnu2DPlaneStrain #undef __FUNCT__ #define __FUNCT__ "MEF90HookeLawIsoLambdaMu3D" subroutine MEF90HookeLawIsoLambdaMu3D(A, lambda, mu) PetscReal, intent(IN) :: lambda, mu type(Tens4OS3D), intent(OUT) :: A A = 0.0_kr A%XXXX = lambda + mu * 2.0_kr A%XXYY = lambda A%XXZZ = lambda A%XYXY = mu A%XZXZ = mu A%YYYY = lambda + mu * 2.0_kr A%YYZZ = lambda A%YZYZ = mu A%ZZZZ = lambda + mu * 2.0_kr end subroutine MEF90HookeLawIsoLambdaMu3D #undef __FUNCT__ #define __FUNCT__ "MEF90HookeLawIsoENu3D" subroutine MEF90HookeLawIsoENu3D(A, E, nu) PetscReal, intent(IN) :: E, nu type(Tens4OS3D), intent(OUT) :: A real(Kind=Kr) :: Lambda, mu lambda = E * nu / (1.0_kr + nu) / (1 - 2.0_kr * nu) mu = E / (1.0_kr + nu)*.5_kr A = 0.0_kr A%XXXX = lambda + mu * 2.0_kr A%XXYY = lambda A%XXZZ = lambda A%XYXY = mu A%XZXZ = mu A%YYYY = lambda + mu * 2.0_kr A%YYZZ = lambda A%YZYZ = mu A%ZZZZ = lambda + mu * 2.0_kr end subroutine MEF90HookeLawIsoENu3D !! Overloading linear algebra functions with Hookes Laws. #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DSum" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DSum: !!! function MEF90HookesLaw2DSum(A, B) type(MEF90HookesLaw2D), intent(IN) :: A, B type(MEF90HookesLaw2D) :: MEF90HookesLaw2DSum character(len=MEF90MXSTRLEN) :: IOBuffer PetscErrorCode :: ierr if ((A%type == MEF90HookesLawTypeIsotropic) .and. (B%type == MEF90HookesLawTypeIsotropic)) then MEF90HookesLaw2DSum%type = MEF90HookesLawTypeIsotropic if (A%isPlaneStress .eqv. B%isPlaneStress) then MEF90HookesLaw2DSum%lambda = A%lambda + B%lambda MEF90HookesLaw2DSum%mu = A%mu + B%mu MEF90HookesLaw2DSum%isPlaneStress = A%isPlaneStress else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if if (A%isPlaneStress) then MEF90HookesLaw2DSum%PoissonRatio = MEF90HookesLaw2DSum%lambda / (MEF90HookesLaw2DSum%lambda + MEF90HookesLaw2DSum%mu) * 0.5_kr MEF90HookesLaw2DSum%YoungsModulus = 2.0_kr * MEF90HookesLaw2DSum%mu * (1.0_kr + MEF90HookesLaw2DSum%PoissonRatio) MEF90HookesLaw2DSum%BulkModulus = MEF90HookesLaw2DSum%lambda + MEF90HookesLaw2DSum%mu else MEF90HookesLaw2DSum%PoissonRatio = MEF90HookesLaw2DSum%lambda / (MEF90HookesLaw2DSum%lambda + 2.0_kr * MEF90HookesLaw2DSum%mu) * 0.5_kr MEF90HookesLaw2DSum%YoungsModulus = 2.0_kr * MEF90HookesLaw2DSum%mu * (1.0_kr + MEF90HookesLaw2DSum%PoissonRatio) MEF90HookesLaw2DSum%BulkModulus = MEF90HookesLaw2DSum%lambda + MEF90HookesLaw2DSum%mu end if else if ((A%type == MEF90HookesLawTypeFull) .and. (B%type == MEF90HookesLawTypeFull)) then MEF90HookesLaw2DSum%type = MEF90HookesLawTypeFull MEF90HookesLaw2DSum%fullTensor = A%fullTensor + B%fullTensor else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if end function MEF90HookesLaw2DSum #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw3DSum" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! MEF90HookesLaw3DSum: !!! function MEF90HookesLaw3DSum(A, B) type(MEF90HookesLaw3D), intent(IN) :: A, B type(MEF90HookesLaw3D) :: MEF90HookesLaw3DSum character(len=MEF90MXSTRLEN) :: IOBuffer PetscErrorCode :: ierr if ((A%type == MEF90HookesLawTypeIsotropic) .and. (B%type == MEF90HookesLawTypeIsotropic)) then MEF90HookesLaw3DSum%type = MEF90HookesLawTypeIsotropic MEF90HookesLaw3DSum%lambda = A%lambda + B%lambda MEF90HookesLaw3DSum%mu = A%mu + B%mu MEF90HookesLaw3DSum%PoissonRatio = MEF90HookesLaw3DSum%lambda / (MEF90HookesLaw3DSum%lambda + MEF90HookesLaw3DSum%mu) * 0.5_kr MEF90HookesLaw3DSum%YoungsModulus = MEF90HookesLaw3DSum%mu * (3.0_kr * MEF90HookesLaw3DSum%lambda + 2.0_kr * MEF90HookesLaw3DSum%mu) / (MEF90HookesLaw3DSum%lambda + MEF90HookesLaw3DSum%mu) MEF90HookesLaw3DSum%BulkModulus = MEF90HookesLaw3DSum%lambda + MEF90HookesLaw3DSum%mu * 2.0_kr / 3.0_kr else if ((A%type == MEF90HookesLawTypeFull) .and. (B%type == MEF90HookesLawTypeFull)) then MEF90HookesLaw3DSum%type = MEF90HookesLawTypeFull MEF90HookesLaw3DSum%fullTensor = A%fullTensor + B%fullTensor else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if end function MEF90HookesLaw3DSum #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DDiff" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DDiff: !!! function MEF90HookesLaw2DDiff(A, B) type(MEF90HookesLaw2D), intent(IN) :: A, B type(MEF90HookesLaw2D) :: MEF90HookesLaw2DDiff character(len=MEF90MXSTRLEN) :: IOBuffer PetscErrorCode :: ierr if ((A%type == MEF90HookesLawTypeIsotropic) .and. (B%type == MEF90HookesLawTypeIsotropic)) then MEF90HookesLaw2DDiff%type = MEF90HookesLawTypeIsotropic if (A%isPlaneStress .eqv. B%isPlaneStress) then MEF90HookesLaw2DDiff%lambda = A%lambda - B%lambda MEF90HookesLaw2DDiff%mu = A%mu - B%mu MEF90HookesLaw2DDiff%isPlaneStress = A%isPlaneStress else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if if (A%isPlaneStress) then MEF90HookesLaw2DDiff%PoissonRatio = MEF90HookesLaw2DDiff%lambda / (MEF90HookesLaw2DDiff%lambda + MEF90HookesLaw2DDiff%mu) * 0.5_kr MEF90HookesLaw2DDiff%YoungsModulus = 2.0_kr * MEF90HookesLaw2DDiff%mu * (1.0_kr + MEF90HookesLaw2DDiff%PoissonRatio) MEF90HookesLaw2DDiff%BulkModulus = MEF90HookesLaw2DDiff%lambda + MEF90HookesLaw2DDiff%mu else MEF90HookesLaw2DDiff%PoissonRatio = MEF90HookesLaw2DDiff%lambda / (MEF90HookesLaw2DDiff%lambda + 2.0_kr * MEF90HookesLaw2DDiff%mu) * 0.5_kr MEF90HookesLaw2DDiff%YoungsModulus = 2.0_kr * MEF90HookesLaw2DDiff%mu * (1.0_kr + MEF90HookesLaw2DDiff%PoissonRatio) MEF90HookesLaw2DDiff%BulkModulus = MEF90HookesLaw2DDiff%lambda + MEF90HookesLaw2DDiff%mu end if else if ((A%type == MEF90HookesLawTypeFull) .and. (B%type == MEF90HookesLawTypeFull)) then MEF90HookesLaw2DDiff%type = MEF90HookesLawTypeFull MEF90HookesLaw2DDiff%fullTensor = A%fullTensor - B%fullTensor else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if end function MEF90HookesLaw2DDiff #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw3DDiff" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! MEF90HookesLaw3DDiff: !!! function MEF90HookesLaw3DDiff(A, B) type(MEF90HookesLaw3D), intent(IN) :: A, B type(MEF90HookesLaw3D) :: MEF90HookesLaw3DDiff character(len=MEF90MXSTRLEN) :: IOBuffer PetscErrorCode :: ierr if ((A%type == MEF90HookesLawTypeIsotropic) .and. (B%type == MEF90HookesLawTypeIsotropic)) then MEF90HookesLaw3DDiff%type = MEF90HookesLawTypeIsotropic MEF90HookesLaw3DDiff%lambda = A%lambda - B%lambda MEF90HookesLaw3DDiff%mu = A%mu - B%mu MEF90HookesLaw3DDiff%PoissonRatio = MEF90HookesLaw3DDiff%lambda / (MEF90HookesLaw3DDiff%lambda + MEF90HookesLaw3DDiff%mu) * 0.5_kr MEF90HookesLaw3DDiff%YoungsModulus = MEF90HookesLaw3DDiff%mu * (3.0_kr * MEF90HookesLaw3DDiff%lambda + 2.0_kr * MEF90HookesLaw3DDiff%mu) / (MEF90HookesLaw3DDiff%lambda + MEF90HookesLaw3DDiff%mu) MEF90HookesLaw3DDiff%BulkModulus = MEF90HookesLaw3DDiff%lambda + MEF90HookesLaw3DDiff%mu * 2.0_kr / 3.0_kr else if ((A%type == MEF90HookesLawTypeFull) .and. (B%type == MEF90HookesLawTypeFull)) then MEF90HookesLaw3DDiff%type = MEF90HookesLawTypeFull MEF90HookesLaw3DDiff%fullTensor = A%fullTensor - B%fullTensor else write (IOBuffer, *) "Incompatible planar Hooke law type in "//__FUNCT__//'\n' PetscCall(PetscPrintf(PETSC_COMM_SELF, IOBuffer, ierr)) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, IOBuffer) end if end function MEF90HookesLaw3DDiff #undef __FUNCT__ #define __FUNCT__ "ScalarXMEF90HookesLaw2D" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! ScalarXMEF90HookesLaw2D: !!! function ScalarXMEF90HookesLaw2D(t, A) PetscReal, intent(IN) :: t type(MEF90HookesLaw2D), intent(IN) :: A type(MEF90HookesLaw2D) :: ScalarXMEF90HookesLaw2D if (A%type == MEF90HookesLawTypeIsotropic) then ScalarXMEF90HookesLaw2D%type = MEF90HookesLawTypeIsotropic ScalarXMEF90HookesLaw2D%lambda = t * A%lambda ScalarXMEF90HookesLaw2D%mu = t * A%mu ScalarXMEF90HookesLaw2D%isPlaneStress = A%isPlaneStress if (A%isPlaneStress) then ScalarXMEF90HookesLaw2D%PoissonRatio = ScalarXMEF90HookesLaw2D%lambda / (ScalarXMEF90HookesLaw2D%lambda + ScalarXMEF90HookesLaw2D%mu) * 0.5_kr ScalarXMEF90HookesLaw2D%YoungsModulus = 2.0_kr * ScalarXMEF90HookesLaw2D%mu * (1.0_kr + ScalarXMEF90HookesLaw2D%PoissonRatio) ScalarXMEF90HookesLaw2D%BulkModulus = ScalarXMEF90HookesLaw2D%lambda + ScalarXMEF90HookesLaw2D%mu else ScalarXMEF90HookesLaw2D%PoissonRatio = ScalarXMEF90HookesLaw2D%lambda / (ScalarXMEF90HookesLaw2D%lambda + 2.0_kr * ScalarXMEF90HookesLaw2D%mu) * 0.5_kr ScalarXMEF90HookesLaw2D%YoungsModulus = 2.0_kr * ScalarXMEF90HookesLaw2D%mu * (1.0_kr + ScalarXMEF90HookesLaw2D%PoissonRatio) ScalarXMEF90HookesLaw2D%BulkModulus = ScalarXMEF90HookesLaw2D%lambda + ScalarXMEF90HookesLaw2D%mu end if else ScalarXMEF90HookesLaw2D%type = MEF90HookesLawTypeFull ScalarXMEF90HookesLaw2D%fullTensor = t * A%fullTensor end if end function ScalarXMEF90HookesLaw2D #undef __FUNCT__ #define __FUNCT__ "ScalarXMEF90HookesLaw3D" !!! author: Blaise Bourdin (2020, bourdin@lsu.edu) !!! !!! ScalarXMEF90HookesLaw3D: !!! function ScalarXMEF90HookesLaw3D(t, A) PetscReal, intent(IN) :: t type(MEF90HookesLaw3D), intent(IN) :: A type(MEF90HookesLaw3D) :: ScalarXMEF90HookesLaw3D if (A%type == MEF90HookesLawTypeIsotropic) then ScalarXMEF90HookesLaw3D%type = MEF90HookesLawTypeIsotropic ScalarXMEF90HookesLaw3D%lambda = t * A%lambda ScalarXMEF90HookesLaw3D%mu = t * A%mu ScalarXMEF90HookesLaw3D%PoissonRatio = ScalarXMEF90HookesLaw3D%lambda / (ScalarXMEF90HookesLaw3D%lambda + ScalarXMEF90HookesLaw3D%mu) * 0.5_kr ScalarXMEF90HookesLaw3D%YoungsModulus = ScalarXMEF90HookesLaw3D%mu * (3.0_kr * ScalarXMEF90HookesLaw3D%lambda + 2.0_kr * ScalarXMEF90HookesLaw3D%mu) / (ScalarXMEF90HookesLaw3D%lambda + ScalarXMEF90HookesLaw3D%mu) ScalarXMEF90HookesLaw3D%BulkModulus = ScalarXMEF90HookesLaw3D%lambda + ScalarXMEF90HookesLaw3D%mu * 2.0_kr / 3.0_kr else ScalarXMEF90HookesLaw3D%type = MEF90HookesLawTypeFull ScalarXMEF90HookesLaw3D%fullTensor = t * A%fullTensor end if end function ScalarXMEF90HookesLaw3D #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DXMatS2D" !!! author: Blaise Bourdin (2016, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DXMatS2D: !!! function MEF90HookesLaw2DXMatS2D(A, X) type(MEF90HookesLaw2D), intent(IN) :: A type(MatS2D), intent(IN) :: X type(MatS2D) :: MEF90HookesLaw2DXMatS2D real(Kind=Kr) :: C1, C2 select case (A%type) case (MEF90HookesLawTypeIsotropic) C1 = A%lambda + 2.0_kr * A%mu C2 = 2.0_kr * A%mu MEF90HookesLaw2DXMatS2D%XX = C1 * X%XX + A%lambda * X%YY MEF90HookesLaw2DXMatS2D%YY = A%lambda * X%XX + C1 * X%YY MEF90HookesLaw2DXMatS2D%XY = C2 * X%XY case (MEF90HookesLawTypeFull) MEF90HookesLaw2DXMatS2D = A%fullTensor * X end select end function MEF90HookesLaw2DXMatS2D #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DXMatS3D" !!! author: Blaise Bourdin (2016, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DXMatS3D: !!! function MEF90HookesLaw3DXMatS3D(A, X) type(MEF90HookesLaw3D), intent(IN) :: A type(MatS3D), intent(IN) :: X type(MatS3D) :: MEF90HookesLaw3DXMatS3D real(Kind=Kr) :: C1, C2 select case (A%type) case (MEF90HookesLawTypeIsotropic) C1 = A%lambda + 2.0_kr * A%mu C2 = 2.0_kr * A%mu MEF90HookesLaw3DXMatS3D%XX = C1 * X%XX + A%lambda * X%YY + A%lambda * X%ZZ MEF90HookesLaw3DXMatS3D%YY = A%lambda * X%XX + C1 * X%YY + A%lambda * X%ZZ MEF90HookesLaw3DXMatS3D%ZZ = A%lambda * X%XX + A%lambda * X%YY + C1 * X%ZZ MEF90HookesLaw3DXMatS3D%YZ = C2 * X%YZ MEF90HookesLaw3DXMatS3D%XZ = C2 * X%XZ MEF90HookesLaw3DXMatS3D%XY = C2 * X%XY case (MEF90HookesLawTypeFull) MEF90HookesLaw3DXMatS3D = A%fullTensor * X end select end function MEF90HookesLaw3DXMatS3D #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DXMat2D" !!! author: Blaise Bourdin (2016, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DXMat2D: !!! function MEF90HookesLaw2DXMat2D(A, X) type(MEF90HookesLaw2D), intent(IN) :: A type(Mat2D), intent(IN) :: X type(Mat2D) :: MEF90HookesLaw2DXMat2D real(Kind=Kr) :: C1, C2 select case (A%type) case (MEF90HookesLawTypeIsotropic) C1 = A%lambda + 2.0_kr * A%mu C2 = 2.0_kr * A%mu MEF90HookesLaw2DXMat2D%XX = C1 * X%XX + A%lambda * X%YY MEF90HookesLaw2DXMat2D%XY = C2 * X%XY MEF90HookesLaw2DXMat2D%YY = A%lambda * X%XX + C1 * X%YY MEF90HookesLaw2DXMat2D%YX = C2 * X%YX case (MEF90HookesLawTypeFull) MEF90HookesLaw2DXMat2D = A%fullTensor * X end select end function MEF90HookesLaw2DXMat2D #undef __FUNCT__ #define __FUNCT__ "MEF90HookesLaw2DXMat3D" !!! author: Blaise Bourdin (2016, bourdin@lsu.edu) !!! !!! MEF90HookesLaw2DXMat3D: !!! function MEF90HookesLaw3DXMat3D(A, X) type(MEF90HookesLaw3D), intent(IN) :: A type(Mat3D), intent(IN) :: X type(Mat3D) :: MEF90HookesLaw3DXMat3D real(Kind=Kr) :: C1, C2 select case (A%type) case (MEF90HookesLawTypeIsotropic) C1 = A%lambda + 2.0_kr * A%mu C2 = 2.0_kr * A%mu MEF90HookesLaw3DXMat3D%XX = C1 * X%XX + A%lambda * X%YY + A%lambda * X%ZZ MEF90HookesLaw3DXMat3D%XY = C2 * X%XY MEF90HookesLaw3DXMat3D%XZ = C2 * X%XZ MEF90HookesLaw3DXMat3D%YX = C2 * X%YX MEF90HookesLaw3DXMat3D%YY = A%lambda * X%XX + C1 * X%YY + A%lambda * X%ZZ MEF90HookesLaw3DXMat3D%YZ = C2 * X%YZ MEF90HookesLaw3DXMat3D%ZX = C2 * X%ZX MEF90HookesLaw3DXMat3D%ZY = C2 * X%ZY MEF90HookesLaw3DXMat3D%ZZ = A%lambda * X%XX + A%lambda * X%YY + C1 * X%ZZ case (MEF90HookesLawTypeFull) MEF90HookesLaw3DXMat3D = A%fullTensor * X end select end function MEF90HookesLaw3DXMat3D end module m_MEF90_Materials