module m_MEF90_LinAlg_class
   implicit none(type, external)
   private
   public :: mef90Vect
   public :: mef90Mat
   public :: mef90Tens4OS

!!! author: Blaise Bourdin (2020, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  m_MEF90_LinAlg_class: Starting a OO implementation of the basic classes in m_MEF90_LinAlg
!!!

   type :: mef90Vect
   end type mef90Vect

   type :: mef90Mat
   end type mef90Mat

   type :: mef90Tens4OS
   end type mef90Tens4OS
end module m_MEF90_LinAlg_class

module m_MEF90_LinAlg
#include "petsc/finclude/petsc.h"
   use m_MEF90_LinAlg_class
   use m_MEF90_Parameters
   use m_MEF90_Utils

   implicit none(type)

   external :: dsyevd
   external :: dgetrf
   external :: dgetri

   type, extends(mef90Vect) :: Vect2D
      PetscReal          :: X = 0.0_Kr
      PetscReal          :: Y = 0.0_Kr
   end type Vect2D

   type, extends(mef90Vect) :: Vect3D
      PetscReal          :: X = 0.0_Kr
      PetscReal          :: Y = 0.0_Kr
      PetscReal          :: Z = 0.0_Kr
   end type Vect3D

   type, extends(mef90Mat) :: Mat2D
      PetscReal          :: XX = 0.0_Kr
      PetscReal          :: XY = 0.0_Kr
      PetscReal          :: YX = 0.0_Kr
      PetscReal          :: YY = 0.0_Kr
   end type Mat2D

   type, extends(mef90Mat) :: MatS2D
      PetscReal          :: XX = 0.0_Kr
      PetscReal          :: YY = 0.0_Kr
      PetscReal          :: XY = 0.0_Kr
   end type MatS2D

   type, extends(mef90Mat) :: Mat3D
      PetscReal          :: XX = 0.0_Kr
      PetscReal          :: XY = 0.0_Kr
      PetscReal          :: XZ = 0.0_Kr
      PetscReal          :: YX = 0.0_Kr
      PetscReal          :: YY = 0.0_Kr
      PetscReal          :: YZ = 0.0_Kr
      PetscReal          :: ZX = 0.0_Kr
      PetscReal          :: ZY = 0.0_Kr
      PetscReal          :: ZZ = 0.0_Kr
   end type Mat3D

   type, extends(mef90Mat) :: MatS3D
      PetscReal          :: XX = 0.0_Kr
      PetscReal          :: YY = 0.0_Kr
      PetscReal          :: ZZ = 0.0_Kr
      PetscReal          :: YZ = 0.0_Kr
      PetscReal          :: XZ = 0.0_Kr
      PetscReal          :: XY = 0.0_Kr
   end type MatS3D

 !! After much hesitation,
 !! - the terms are numbered in alphabetical order (i.e. XXYX and not XYXX)
 !! - the terms are stored in alphabetical order
 !! 2014-07: Changed ordering to rows of the upper triangular part
 !!          and naming to be consistent with Voigt notations
   type, extends(mef90Tens4OS) :: Tens4OS2D
      PetscReal          :: XXXX, XXYY, XXXY
      PetscReal          ::       YYYY, YYXY
      PetscReal          ::             XYXY

   end type Tens4OS2D

   type, extends(mef90Tens4OS) :: Tens4OS3D
      PetscReal          :: XXXX, XXYY, XXZZ, XXYZ, XXXZ, XXXY
      PetscReal          ::       YYYY, YYZZ, YYYZ, YYXZ, YYXY
      PetscReal          ::             ZZZZ, ZZYZ, ZZXZ, ZZXY
      PetscReal          ::                   YZYZ, YZXZ, YZXY
      PetscReal          ::                         XZXZ, XZXY
      PetscReal          ::                               XYXY
   end type Tens4OS3D

   type(Vect2D), parameter    :: MEF90Vect2De1 = Vect2D(1.0_kr, 0.0_kr)
   type(Vect2D), parameter    :: MEF90Vect2De2 = Vect2D(0.0_kr, 1.0_kr)

   type(Vect3D), parameter    :: MEF90Vect3De1 = Vect3D(1.0_kr, 0.0_kr, 0.0_kr)
   type(Vect3D), parameter    :: MEF90Vect3De2 = Vect3D(0.0_kr, 1.0_kr, 0.0_kr)
   type(Vect3D), parameter    :: MEF90Vect3De3 = Vect3D(0.0_kr, 0.0_kr, 1.0_kr)

   type(Mat2D), parameter     :: MEF90Mat2DIdentity = Mat2D(1.0_kr, 0.0_kr, &
                                                            0.0_kr, 1.0_kr)
   type(MatS2D), parameter    :: MEF90MatS2DIdentity = MatS2D(1.0_kr, 1.0_kr, 0.0_kr)
   type(Mat3D), parameter     :: MEF90Mat3DIdentity = Mat3D(1.0_kr, 0.0_kr, 0.0_kr, &
                                                            0.0_kr, 1.0_kr, 0.0_kr, &
                                                            0.0_kr, 0.0_kr, 1.0_kr)
   type(MatS3D), parameter    :: MEF90MatS3DIdentity = MatS3D(1.0_kr, 1.0_kr, 1.0_kr, &
                                                              0.0_kr, 0.0_kr, 0.0_kr)

   type(Tens4OS2D), parameter :: MEF90Tens4OS2DIdentity = Tens4OS2D(1.0_kr, 0.0_kr, 0.0_kr, &
                                                                    1.0_kr, 0.0_kr, &
                                                                    0.5_kr)
   type(Tens4OS3D), parameter :: MEF90Tens4OS3DIdentity = Tens4OS3D(1.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, &
                                                                    1.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, &
                                                                    1.0_kr, 0.0_kr, 0.0_kr, 0.0_kr, &
                                                                    0.5_kr, 0.0_kr, 0.0_kr, &
                                                                    0.5_kr, 0.0_kr, &
                                                                    0.5_kr)

   interface operator(+)
      module procedure SumVect2D, SumVect3D, SumMat2D, SumMat3D, SumMatS2D, SumMatS3D, SumTens4OS2D, SumTens4OS3D
   end interface

   interface operator(-)
      module procedure DifVect2D, DifVect3D, DifMat2D, DifMat3D, DifMatS2D, DifMatS3D, DifTens4OS2D, DifTens4OS3D
   end interface

   interface operator(*)
      module procedure DbleXVect2D, Vect2DXDble, DbleXVect3D, Vect3DXDble, &
         DbleXMat2D, Mat2DXDble, DbleXMat3D, Mat3DXDble, &
         DbleXMatS2D, MatS2DXDble, DbleXMatS3D, MatS3DXDble, &
         MatXVect2D, MatXVect3D, MatXVect2DS, MatXVect3DS, &
         DbleXTens4OS2D, Tens4OS2DXDble, Tens4OS2DXMatS2D, Tens4OS2DXMat2D, &
         DbleXTens4OS3D, Tens4OS3DXDble, Tens4OS3DXMatS3D, Tens4OS3DXMat3D, &
         Mat2DXMat2D, Mat2DXMatS2D, MatS2DXMatS2D, Mat3DXMat3D, Mat3DXMatS3D, MatS3DXMatS3D, &
         DotP2D, DotP3D
   end interface

   interface operator(/)
      module procedure Vect2DQuot, Vect3DQuot, Mat2DQuot, Mat3DQuot, MatS2DQuot, MatS3DQuot, Tens4OS2DQuot, Tens4OS3DQuot
   end interface

   interface operator(.DotP.)
      module procedure DotP2D, DotP3D, ContP2D, ContP3D, ContP2DS, ContP3DS, &
         Mat2DDotMatS2D, MatS2DDotMat2D, Mat3DDotMatS3D, MatS3DDotMat3D
   end interface

   interface operator(.CrossP.)
      module procedure CrossP3D
   end interface

   interface Transpose
      module procedure Transpose2D, Transpose3D
   end interface

   interface Invert
      module procedure InvertMat2D, InvertMatS2D, InvertMat3D, InvertMatS3D, InvertTens4OS2D, InvertTens4OS3D
   end interface

   interface operator(.TensP.)
      module procedure TensPVect2D, TensPVect3D
   end interface

   interface operator(.SymP.)
      module procedure SymPVect2D, SymPVect3D, SymPMatS2D, SymPMatS3D
   end interface

   interface operator(.oDot.)
      module procedure oDotMatS2D, oDotMatS3D
   end interface

   interface Trace
      module procedure Trace2D, Trace3D, Trace2DS, Trace3DS
   end interface

   interface Det
      module procedure DetMat2D, DetMatS2D, DetMat3D, DetMatS3D
   end interface

   interface assignment(=)
      module procedure Vect2D_Get_Real, Vect3D_Get_Real, &
         Vect2D_Get_VectR, Vect3D_Get_VectR, &
         Vect2DEQ, Vect3DEQ, Mat2D_Get_Real, Mat3D_Get_Real, &
         Mat2DEQ, Mat3DEQ, MatS2D_Get_Real, MatS3D_Get_Real, &
         MatS2DEQ, MatS3DEQ, MatS2D_Get_VectR, MatS3D_Get_VectR, &
         VectR_Get_MatS2D, VectR_Get_MatS3D, &
         Tens4OS2D_Get_Real, Tens4OS3D_Get_Real, &
         Tens4OS2DToArray, ArrayToTens4OS2D, &
         Tens4OS3DToArray, ArrayToTens4OS3D, &
         Mat2DGetArray, Mat3DGetArray, &
         MatS2DGetArray, MatS3DGetArray, &
         ArrayGetMat2D, ArrayGetMat3D, &
         ArrayGetMatS2D, ArrayGetMatS3D, &
         MatS3DToMat3D, MatS2DToMat2D, &
         Mat3DToMatS3D, Mat2DToMatS2D, &
         Tens4OS2D2Array4, Tens4OS3D2Array4, &
         Array42Tens4OS2D, Array42Tens4OS3D
   end interface

   interface Symmetrize
      module procedure Symmetrize2D, Symmetrize3D
   end interface

   interface DeviatoricPart
      module procedure DeviatoricPart2D, DeviatoricPart2DS, DeviatoricPart3D, DeviatoricPart3DS
   end interface

   interface HydrostaticPart
      module procedure HydrostaticPart2D, HydrostaticPart2DS, HydrostaticPart3D, HydrostaticPart3DS
   end interface

   interface MEF90MatRaRt
      module procedure RaRtMat2D, RaRtMatS2D, RaRtMat3D, RaRtMatS3D
   end interface

   interface MEF90MatRtaR
      module procedure RtaRMat2D, RtaRMatS2D, RtaRMat3D, RtaRMatS3D
   end interface

   interface Norm
      module procedure Vect2DNorm, Vect3DNorm, Mat2DNorm, MatS2DNorm, Mat3DNorm, MatS3DNorm
      ! Tens4OS2DNorm,Tens4OS3DNorm
   end interface

   interface simplexNormal
      module procedure simplexNormal2D, simplexNormal3D
   end interface simplexNormal

   interface sqrt
      module procedure Tens4OS2DSquareRoot, Tens4OS3DSquareRoot
   end interface

   interface SpectralDecomposition
      module procedure MatS2DSpectralDecomposition, MatS3DSpectralDecomposition
   end interface

   interface Diagonalize
      !!! Diagonalize(A,P,D) returns P,D such that A = P D P^{-1}
      !!! The diagonal entries of D are sorted in increasing order.
      module procedure MatS2DEigenVectorValues, MatS3DEigenVectorValues
   end interface

   interface Moment
      module procedure Mat2DMoment, MatS2DMoment, Mat3DMoment, MatS3DMoment
   end interface

   interface Tens4OSTransform
      module procedure Tens4OS2DTransform, Tens4OS3DTransform
   end interface

   !Interface MatSymToMat
   !   Module Procedure MatS3DToMat3D,MatS2DToMat2D
   !End Interface

!!$  Type(Vect2D),Parameter       :: e1_2D = [1.0_Kr,0.0_Kr]
!!$  Type(Vect2D),Parameter       :: e2_2D = [0.0_Kr,1.0_Kr]
!!$
!!$  Type(Vect3D),Parameter       :: e1_3D = [1.0_Kr,0.0_Kr,0.0_Kr]
!!$  Type(Vect3D),Parameter       :: e2_3D = [0.0_Kr,1.0_Kr,0.0_Kr]
!!$  Type(Vect3D),Parameter       :: e3_3D = [0.0_Kr,0.0_Kr,1.0_Kr]

contains
   function SumVect2D(V1, V2)
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D), intent(IN)                     :: V2
      type(Vect2D)                                :: SumVect2D

      SumVect2D%X = V1%X + V2%X
      SumVect2D%Y = V1%Y + V2%Y
   end function SumVect2D

   function SumVect3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      type(Vect3D)                                :: SumVect3D

      SumVect3D%X = V1%X + V2%X
      SumVect3D%Y = V1%Y + V2%Y
      SumVect3D%Z = V1%Z + V2%Z
   end function SumVect3D

   function SumMat2D(M1, M2)
      type(Mat2D), intent(IN)                      :: M1, M2
      type(Mat2D)                                 :: SumMat2D

      SumMat2D%XX = M1%XX + M2%XX
      SumMat2D%XY = M1%XY + M2%XY
      SumMat2D%YX = M1%YX + M2%YX
      SumMat2D%YY = M1%YY + M2%YY
   end function SumMat2D

   function SumMatS2D(M1, M2)
      type(MatS2D), intent(IN)                     :: M1, M2
      type(MatS2D)                                :: SumMatS2D

      SumMatS2D%XX = M1%XX + M2%XX
      SumMatS2D%XY = M1%XY + M2%XY
      SumMatS2D%YY = M1%YY + M2%YY
   end function SumMatS2D

   function SumMat3D(M1, M2)
      type(Mat3D), intent(IN)                      :: M1, M2
      type(Mat3D)                                 :: SumMat3D

      SumMat3D%XX = M1%XX + M2%XX
      SumMat3D%XY = M1%XY + M2%XY
      SumMat3D%XZ = M1%XZ + M2%XZ
      SumMat3D%YX = M1%YX + M2%YX
      SumMat3D%YY = M1%YY + M2%YY
      SumMat3D%YZ = M1%YZ + M2%YZ
      SumMat3D%ZX = M1%ZX + M2%ZX
      SumMat3D%ZY = M1%ZY + M2%ZY
      SumMat3D%ZZ = M1%ZZ + M2%ZZ
   end function SumMat3D

   function SumMatS3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1, M2
      type(MatS3D)                                :: SumMatS3D

      SumMatS3D%XX = M1%XX + M2%XX
      SumMatS3D%YY = M1%YY + M2%YY
      SumMatS3D%ZZ = M1%ZZ + M2%ZZ
      SumMatS3D%YZ = M1%YZ + M2%YZ
      SumMatS3D%XZ = M1%XZ + M2%XZ
      SumMatS3D%XY = M1%XY + M2%XY
   end function SumMatS3D

   function SumTens4OS2D(T1, T2)
      type(Tens4OS2D), intent(IN)                  :: T1, T2
      type(Tens4OS2D)                             :: SumTens4OS2D

      SumTens4OS2D%XXXX = T1%XXXX + T2%XXXX
      SumTens4OS2D%XXYY = T1%XXYY + T2%XXYY
      SumTens4OS2D%XXXY = T1%XXXY + T2%XXXY

      SumTens4OS2D%YYYY = T1%YYYY + T2%YYYY
      SumTens4OS2D%YYXY = T1%YYXY + T2%YYXY

      SumTens4OS2D%XYXY = T1%XYXY + T2%XYXY

   end function SumTens4OS2D

   function SumTens4OS3D(T1, T2)
      type(Tens4OS3D), intent(IN)                  :: T1, T2
      type(Tens4OS3D)                             :: SumTens4OS3D

      SumTens4OS3D%XXXX = T1%XXXX + T2%XXXX
      SumTens4OS3D%XXYY = T1%XXYY + T2%XXYY
      SumTens4OS3D%XXZZ = T1%XXZZ + T2%XXZZ
      SumTens4OS3D%XXYZ = T1%XXYZ + T2%XXYZ
      SumTens4OS3D%XXXZ = T1%XXXZ + T2%XXXZ
      SumTens4OS3D%XXXY = T1%XXXY + T2%XXXY

      SumTens4OS3D%YYYY = T1%YYYY + T2%YYYY
      SumTens4OS3D%YYZZ = T1%YYZZ + T2%YYZZ
      SumTens4OS3D%YYYZ = T1%YYYZ + T2%YYYZ
      SumTens4OS3D%YYXZ = T1%YYXZ + T2%YYXZ
      SumTens4OS3D%YYXY = T1%YYXY + T2%YYXY

      SumTens4OS3D%ZZZZ = T1%ZZZZ + T2%ZZZZ
      SumTens4OS3D%ZZYZ = T1%ZZYZ + T2%ZZYZ
      SumTens4OS3D%ZZXZ = T1%ZZXZ + T2%ZZXZ
      SumTens4OS3D%ZZXY = T1%ZZXY + T2%ZZXY

      SumTens4OS3D%YZYZ = T1%YZYZ + T2%YZYZ
      SumTens4OS3D%YZXZ = T1%YZXZ + T2%YZXZ
      SumTens4OS3D%YZXY = T1%YZXY + T2%YZXY

      SumTens4OS3D%XZXZ = T1%XZXZ + T2%XZXZ
      SumTens4OS3D%XZXY = T1%XZXY + T2%XZXY

      SumTens4OS3D%XYXY = T1%XYXY + T2%XYXY
   end function SumTens4OS3D

   ! Overloading "-"
   function DifVect2D(V1, V2)
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D), intent(IN)                     :: V2
      type(Vect2D)                                 :: DifVect2D

      DifVect2D%X = V1%X - V2%X
      DifVect2D%Y = V1%Y - V2%Y
   end function DifVect2D

   function DifVect3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      type(Vect3D)                                 :: DifVect3D

      DifVect3D%X = V1%X - V2%X
      DifVect3D%Y = V1%Y - V2%Y
      DifVect3D%Z = V1%Z - V2%Z
   end function DifVect3D

   function DifMat2D(M1, M2)
      type(Mat2D), intent(IN)                      :: M1, M2
      type(Mat2D)                                  :: DifMat2D

      DifMat2D%XX = M1%XX - M2%XX
      DifMat2D%XY = M1%XY - M2%XY
      DifMat2D%YX = M1%YX - M2%YX
      DifMat2D%YY = M1%YY - M2%YY
   end function DifMat2D

   function DifMatS2D(M1, M2)
      type(MatS2D), intent(IN)                      :: M1, M2
      type(MatS2D)                                  :: DifMatS2D

      DifMatS2D%XX = M1%XX - M2%XX
      DifMatS2D%YY = M1%YY - M2%YY
      DifMatS2D%XY = M1%XY - M2%XY
   end function DifMatS2D

   function DifMat3D(M1, M2)
      type(Mat3D), intent(IN)                      :: M1, M2
      type(Mat3D)                                  :: DifMat3D

      DifMat3D%XX = M1%XX - M2%XX
      DifMat3D%XY = M1%XY - M2%XY
      DifMat3D%XZ = M1%XZ - M2%XZ
      DifMat3D%YX = M1%YX - M2%YX
      DifMat3D%YY = M1%YY - M2%YY
      DifMat3D%YZ = M1%YZ - M2%YZ
      DifMat3D%ZX = M1%ZX - M2%ZX
      DifMat3D%ZY = M1%ZY - M2%ZY
      DifMat3D%ZZ = M1%ZZ - M2%ZZ
   end function DifMat3D

   function DifMatS3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1, M2
      type(MatS3D)                                 :: DifMatS3D

      DifMatS3D%XX = M1%XX - M2%XX
      DifMatS3D%YY = M1%YY - M2%YY
      DifMatS3D%ZZ = M1%ZZ - M2%ZZ
      DifMatS3D%YZ = M1%YZ - M2%YZ
      DifMatS3D%XZ = M1%XZ - M2%XZ
      DifMatS3D%XY = M1%XY - M2%XY
   end function DifMatS3D

   function DifTens4OS2D(T1, T2)
      type(Tens4OS2D), intent(IN)                  :: T1, T2
      type(Tens4OS2D)                              :: DifTens4OS2D

      DifTens4OS2D%XXXX = T1%XXXX - T2%XXXX
      DifTens4OS2D%XXYY = T1%XXYY - T2%XXYY
      DifTens4OS2D%XXXY = T1%XXXY - T2%XXXY

      DifTens4OS2D%YYYY = T1%YYYY - T2%YYYY
      DifTens4OS2D%YYXY = T1%YYXY - T2%YYXY

      DifTens4OS2D%XYXY = T1%XYXY - T2%XYXY
   end function DifTens4OS2D

   function DifTens4OS3D(T1, T2)
      type(Tens4OS3D), intent(IN)                  :: T1, T2
      type(Tens4OS3D)                              :: DifTens4OS3D

      DifTens4OS3D%XXXX = T1%XXXX - T2%XXXX
      DifTens4OS3D%XXYY = T1%XXYY - T2%XXYY
      DifTens4OS3D%XXZZ = T1%XXZZ - T2%XXZZ
      DifTens4OS3D%XXYZ = T1%XXYZ - T2%XXYZ
      DifTens4OS3D%XXXZ = T1%XXXZ - T2%XXXZ
      DifTens4OS3D%XXXY = T1%XXXY - T2%XXXY

      DifTens4OS3D%YYYY = T1%YYYY - T2%YYYY
      DifTens4OS3D%YYZZ = T1%YYZZ - T2%YYZZ
      DifTens4OS3D%YYYZ = T1%YYYZ - T2%YYYZ
      DifTens4OS3D%YYXZ = T1%YYXZ - T2%YYXZ
      DifTens4OS3D%YYXY = T1%YYXY - T2%YYXY

      DifTens4OS3D%ZZZZ = T1%ZZZZ - T2%ZZZZ
      DifTens4OS3D%ZZYZ = T1%ZZYZ - T2%ZZYZ
      DifTens4OS3D%ZZXZ = T1%ZZXZ - T2%ZZXZ
      DifTens4OS3D%ZZXY = T1%ZZXY - T2%ZZXY

      DifTens4OS3D%YZYZ = T1%YZYZ - T2%YZYZ
      DifTens4OS3D%YZXZ = T1%YZXZ - T2%YZXZ
      DifTens4OS3D%YZXY = T1%YZXY - T2%YZXY

      DifTens4OS3D%XZXZ = T1%XZXZ - T2%XZXZ
      DifTens4OS3D%XZXY = T1%XZXY - T2%XZXY

      DifTens4OS3D%XYXY = T1%XYXY - T2%XYXY
   end function DifTens4OS3D

   ! Overloading "*"
   function DbleXVect2D(D1, V1)
      PetscReal, intent(IN)                        :: D1
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D)                                 :: DbleXVect2D

      DbleXVect2D%X = D1 * V1%X
      DbleXVect2D%Y = D1 * V1%Y
   end function DbleXVect2D

   function Vect2DXDble(V1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D)                                 :: Vect2DXDble

      Vect2DXDble%X = D1 * V1%X
      Vect2DXDble%Y = D1 * V1%Y
   end function Vect2DXDble

   function DbleXVect3D(D1, V1)
      PetscReal, intent(IN)                        :: D1
      type(Vect3D), intent(IN)                     :: V1
      type(Vect3D)                                 :: DbleXVect3D

      DbleXVect3D%X = D1 * V1%X
      DbleXVect3D%Y = D1 * V1%Y
      DbleXVect3D%Z = D1 * V1%Z
   end function DbleXVect3D

   function Vect3DXDble(V1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Vect3D), intent(IN)                     :: V1
      type(Vect3D)                                 :: Vect3DXDble

      Vect3DXDble%X = D1 * V1%X
      Vect3DXDble%Y = D1 * V1%Y
      Vect3DXDble%Z = D1 * V1%Z
   end function Vect3DXDble

   function DbleXMat2D(D1, M1)
      PetscReal, intent(IN)                        :: D1
      type(Mat2D), intent(IN)                      :: M1
      type(Mat2D)                                  :: DbleXMat2D

      DbleXMat2D%XX = D1 * M1%XX
      DbleXMat2D%XY = D1 * M1%XY
      DbleXMat2D%YX = D1 * M1%YX
      DbleXMat2D%YY = D1 * M1%YY
   end function DbleXMat2D

   function Mat2DXDble(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Mat2D), intent(IN)                      :: M1
      type(Mat2D)                                  :: Mat2DXDble

      Mat2DXDble%XX = D1 * M1%XX
      Mat2DXDble%XY = D1 * M1%XY
      Mat2DXDble%YX = D1 * M1%YX
      Mat2DXDble%YY = D1 * M1%YY
   end function Mat2DXDble

   function DbleXMatS2D(D1, M1)
      PetscReal, intent(IN)                        :: D1
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: DbleXMatS2D

      DbleXMatS2D%XX = D1 * M1%XX
      DbleXMatS2D%YY = D1 * M1%YY
      DbleXMatS2D%XY = D1 * M1%XY
   end function DbleXMatS2D

   function MatS2DXDble(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: MatS2DXDble

      MatS2DXDble%XX = D1 * M1%XX
      MatS2DXDble%YY = D1 * M1%YY
      MatS2DXDble%XY = D1 * M1%XY
   end function MatS2DXDble

   function DbleXMat3D(D1, M1)
      PetscReal, intent(IN)                        :: D1
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: DbleXMat3D

      DbleXMat3D%XX = D1 * M1%XX
      DbleXMat3D%XY = D1 * M1%XY
      DbleXMat3D%XZ = D1 * M1%XZ
      DbleXMat3D%YX = D1 * M1%YX
      DbleXMat3D%YY = D1 * M1%YY
      DbleXMat3D%YZ = D1 * M1%YZ
      DbleXMat3D%ZX = D1 * M1%ZX
      DbleXMat3D%ZY = D1 * M1%ZY
      DbleXMat3D%ZZ = D1 * M1%ZZ
   end function DbleXMat3D

   function Mat3DXDble(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: Mat3DXDble

      Mat3DXDble%XX = D1 * M1%XX
      Mat3DXDble%XY = D1 * M1%XY
      Mat3DXDble%XZ = D1 * M1%XZ
      Mat3DXDble%YX = D1 * M1%YX
      Mat3DXDble%YY = D1 * M1%YY
      Mat3DXDble%YZ = D1 * M1%YZ
      Mat3DXDble%ZX = D1 * M1%ZX
      Mat3DXDble%ZY = D1 * M1%ZY
      Mat3DXDble%ZZ = D1 * M1%ZZ
   end function Mat3DXDble

   function DbleXMatS3D(D1, M1)
      PetscReal, intent(IN)                        :: D1
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: DbleXMatS3D

      DbleXMatS3D%XX = D1 * M1%XX
      DbleXMatS3D%YY = D1 * M1%YY
      DbleXMatS3D%ZZ = D1 * M1%ZZ
      DbleXMatS3D%YZ = D1 * M1%YZ
      DbleXMatS3D%XZ = D1 * M1%XZ
      DbleXMatS3D%XY = D1 * M1%XY
   end function DbleXMatS3D

   function MatS3DXDble(M1, D1)
      real(Kind=Kr), intent(IN)                    :: D1
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: MatS3DXDble

      MatS3DXDble%XX = D1 * M1%XX
      MatS3DXDble%YY = D1 * M1%YY
      MatS3DXDble%ZZ = D1 * M1%ZZ
      MatS3DXDble%YZ = D1 * M1%YZ
      MatS3DXDble%XZ = D1 * M1%XZ
      MatS3DXDble%XY = D1 * M1%XY
   end function MatS3DXDble

   function DbleXTens4OS2D(D1, T1)
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS2D), intent(IN)                  :: T1
      type(Tens4OS2D)                              :: DbleXTens4OS2D

      DbleXTens4OS2D%XXXX = D1 * T1%XXXX
      DbleXTens4OS2D%XXXY = D1 * T1%XXXY
      DbleXTens4OS2D%XXYY = D1 * T1%XXYY

      DbleXTens4OS2D%XYXY = D1 * T1%XYXY
      DbleXTens4OS2D%YYXY = D1 * T1%YYXY

      DbleXTens4OS2D%YYYY = D1 * T1%YYYY
   end function DbleXTens4OS2D

   function DbleXTens4OS3D(D1, T1)
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS3D), intent(IN)                  :: T1
      type(Tens4OS3D)                              :: DbleXTens4OS3D

      DbleXTens4OS3D%XXXX = D1 * T1%XXXX
      DbleXTens4OS3D%XXYY = D1 * T1%XXYY
      DbleXTens4OS3D%XXZZ = D1 * T1%XXZZ
      DbleXTens4OS3D%XXYZ = D1 * T1%XXYZ
      DbleXTens4OS3D%XXXZ = D1 * T1%XXXZ
      DbleXTens4OS3D%XXXY = D1 * T1%XXXY

      DbleXTens4OS3D%YYYY = D1 * T1%YYYY
      DbleXTens4OS3D%YYZZ = D1 * T1%YYZZ
      DbleXTens4OS3D%YYYZ = D1 * T1%YYYZ
      DbleXTens4OS3D%YYXZ = D1 * T1%YYXZ
      DbleXTens4OS3D%YYXY = D1 * T1%YYXY

      DbleXTens4OS3D%ZZZZ = D1 * T1%ZZZZ
      DbleXTens4OS3D%ZZYZ = D1 * T1%ZZYZ
      DbleXTens4OS3D%ZZXZ = D1 * T1%ZZXZ
      DbleXTens4OS3D%ZZXY = D1 * T1%ZZXY

      DbleXTens4OS3D%YZYZ = D1 * T1%YZYZ
      DbleXTens4OS3D%YZXZ = D1 * T1%YZXZ
      DbleXTens4OS3D%YZXY = D1 * T1%YZXY

      DbleXTens4OS3D%XZXZ = D1 * T1%XZXZ
      DbleXTens4OS3D%XZXY = D1 * T1%XZXY

      DbleXTens4OS3D%XYXY = D1 * T1%XYXY
   end function DbleXTens4OS3D

   function Tens4OS2DXDble(T1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS2D), intent(IN)                  :: T1
      type(Tens4OS2D)                              :: Tens4OS2DXDble

      Tens4OS2DXDble%XXXX = D1 * T1%XXXX
      Tens4OS2DXDble%XXXY = D1 * T1%XXXY
      Tens4OS2DXDble%XXYY = D1 * T1%XXYY

      Tens4OS2DXDble%XYXY = D1 * T1%XYXY
      Tens4OS2DXDble%YYXY = D1 * T1%YYXY

      Tens4OS2DXDble%YYYY = D1 * T1%YYYY
   end function Tens4OS2DXDble

   function Tens4OS3DXDble(T1, D1)
      type(Tens4OS3D), intent(IN)                  :: T1
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS3D)                              :: Tens4OS3DXDble

      Tens4OS3DXDble%XXXX = D1 * T1%XXXX
      Tens4OS3DXDble%XXYY = D1 * T1%XXYY
      Tens4OS3DXDble%XXZZ = D1 * T1%XXZZ
      Tens4OS3DXDble%XXYZ = D1 * T1%XXYZ
      Tens4OS3DXDble%XXXZ = D1 * T1%XXXZ
      Tens4OS3DXDble%XXXY = D1 * T1%XXXY

      Tens4OS3DXDble%YYYY = D1 * T1%YYYY
      Tens4OS3DXDble%YYZZ = D1 * T1%YYZZ
      Tens4OS3DXDble%YYYZ = D1 * T1%YYYZ
      Tens4OS3DXDble%YYXZ = D1 * T1%YYXZ
      Tens4OS3DXDble%YYXY = D1 * T1%YYXY

      Tens4OS3DXDble%ZZZZ = D1 * T1%ZZZZ
      Tens4OS3DXDble%ZZYZ = D1 * T1%ZZYZ
      Tens4OS3DXDble%ZZXZ = D1 * T1%ZZXZ
      Tens4OS3DXDble%ZZXY = D1 * T1%ZZXY

      Tens4OS3DXDble%YZYZ = D1 * T1%YZYZ
      Tens4OS3DXDble%YZXZ = D1 * T1%YZXZ
      Tens4OS3DXDble%YZXY = D1 * T1%YZXY

      Tens4OS3DXDble%XZXZ = D1 * T1%XZXZ
      Tens4OS3DXDble%XZXY = D1 * T1%XZXY

      Tens4OS3DXDble%XYXY = D1 * T1%XYXY
   end function Tens4OS3DXDble

   function MatXVect2D(M1, V1)
      type(Mat2D), intent(IN)                      :: M1
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D)                                 :: MatXVect2D

      MatXVect2D%X = M1%XX * V1%X + M1%XY * V1%Y
      MatXVect2D%Y = M1%YX * V1%X + M1%YY * V1%Y
   end function MatXVect2D

   function MatXVect2DS(M1, V1)
      type(MatS2D), intent(IN)                     :: M1
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D)                                 :: MatXVect2DS

      MatXVect2DS%X = M1%XX * V1%X + M1%XY * V1%Y
      MatXVect2DS%Y = M1%XY * V1%X + M1%YY * V1%Y
   end function MatXVect2DS

   function MatXVect3D(M1, V1)
      type(Mat3D), intent(IN)                      :: M1
      type(Vect3D), intent(IN)                     :: V1
      type(Vect3D)                                 :: MatXVect3D

      MatXVect3D%X = M1%XX * V1%X + M1%XY * V1%Y + M1%XZ * V1%Z
      MatXVect3D%Y = M1%YX * V1%X + M1%YY * V1%Y + M1%YZ * V1%Z
      MatXVect3D%Z = M1%ZX * V1%X + M1%ZY * V1%Y + M1%ZZ * V1%Z
   end function MatXVect3D

   function MatXVect3DS(M1, V1)
      type(MatS3D), intent(IN)                     :: M1
      type(Vect3D), intent(IN)                     :: V1
      type(Vect3D)                                 :: MatXVect3DS

      MatXVect3DS%X = M1%XX * V1%X + M1%XY * V1%Y + M1%XZ * V1%Z
      MatXVect3DS%Y = M1%XY * V1%X + M1%YY * V1%Y + M1%YZ * V1%Z
      MatXVect3DS%Z = M1%XZ * V1%X + M1%YZ * V1%Y + M1%ZZ * V1%Z
   end function MatXVect3DS

   function Tens4OS2DXMatS2D(T1, M1)
      type(Tens4OS2D), intent(IN)                  :: T1
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: Tens4OS2DXMatS2D

      Tens4OS2DXMatS2D%XX = T1%XXXX * M1%XX + T1%XXYY * M1%YY + T1%XXXY * M1%XY * 2.0_kr
      Tens4OS2DXMatS2D%YY = T1%XXYY * M1%XX + T1%YYYY * M1%YY + T1%YYXY * M1%XY * 2.0_kr
      Tens4OS2DXMatS2D%XY = T1%XXXY * M1%XX + T1%YYXY * M1%YY + T1%XYXY * M1%XY * 2.0_kr
   end function Tens4OS2DXMatS2D

   function Tens4OS3DXMatS3D(T1, M1)
      type(Tens4OS3D), intent(IN)                  :: T1
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: Tens4OS3DXMatS3D

      Tens4OS3DXMatS3D%XX = T1%XXXX * M1%XX + T1%XXYY * M1%YY + T1%XXZZ * M1%ZZ &
                              + (T1%XXYZ * M1%YZ + T1%XXXZ * M1%XZ + T1%XXXY * M1%XY) * 2.0_kr
      Tens4OS3DXMatS3D%YY = T1%XXYY * M1%XX + T1%YYYY * M1%YY + T1%YYZZ * M1%ZZ &
                              + (T1%YYYZ * M1%YZ + T1%YYXZ * M1%XZ + T1%YYXY * M1%XY) * 2.0_kr
      Tens4OS3DXMatS3D%ZZ = T1%XXZZ * M1%XX + T1%YYZZ * M1%YY + T1%ZZZZ * M1%ZZ &
                              + (T1%ZZYZ * M1%YZ + T1%ZZXZ * M1%XZ + T1%ZZXY * M1%XY) * 2.0_kr
      Tens4OS3DXMatS3D%YZ = T1%XXYZ * M1%XX + T1%YYYZ * M1%YY + T1%ZZYZ * M1%ZZ &
                              + (T1%YZYZ * M1%YZ + T1%YZXZ * M1%XZ + T1%YZXY * M1%XY) * 2.0_kr
      Tens4OS3DXMatS3D%XZ = T1%XXXZ * M1%XX + T1%YYXZ * M1%YY + T1%ZZXZ * M1%ZZ &
                              + (T1%YZXZ * M1%YZ + T1%XZXZ * M1%XZ + T1%XZXY * M1%XY) * 2.0_kr
      Tens4OS3DXMatS3D%XY = T1%XXXY * M1%XX + T1%YYXY * M1%YY + T1%ZZXY * M1%ZZ &
                              + (T1%YZXY * M1%YZ + T1%XZXY * M1%XZ + T1%XYXY * M1%XY) * 2.0_kr
   end function Tens4OS3DXMatS3D

#undef __FUNCT__
#define __FUNCT__ "Tens4OS2DXMat2D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  Tens4OS2DXMat2D:
!!!
   function Tens4OS2DXMat2D(T1, M1)
      type(Tens4OS2D), intent(IN)                  :: T1
      type(Mat2D), intent(IN)                      :: M1
      type(MatS2D)                                 :: Tens4OS2DXMat2D

      Tens4OS2DXMat2D = T1 * symmetrize(M1)
   end function Tens4OS2DXMat2D

#undef __FUNCT__
#define __FUNCT__ "Tens4OS3DXMat3D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  Tens4OS3DXMat3D:
!!!
   function Tens4OS3DXMat3D(T1, M1)
      type(Tens4OS3D), intent(IN)                  :: T1
      type(Mat3D), intent(IN)                      :: M1
      type(MatS3D)                                 :: Tens4OS3DXMat3D

      Tens4OS3DXMat3D = T1 * symmetrize(M1)
   end function Tens4OS3DXMat3D

   function Mat2DXMat2D(M1, M2)
      type(Mat2D), intent(IN)                      :: M1, M2
      type(Mat2D)                                 :: Mat2DXMat2D

      Mat2DXMat2D%XX = M1%XX * M2%XX + M1%XY * M2%YX
      Mat2DXMat2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      Mat2DXMat2D%YX = M1%YX * M2%XX + M1%YY * M2%YX
      Mat2DXMat2D%YY = M1%YX * M2%XY + M1%YY * M2%YY
   end function Mat2DXMat2D

   function Mat2DXMatS2D(M1, M2)
      type(Mat2D), intent(IN)                      :: M1
      type(MatS2D), intent(IN)                     :: M2
      type(Mat2D)                                  :: Mat2DXMatS2D

      Mat2DXMatS2D%XX = M1%XX * M2%XX + M1%XY * M2%XY
      Mat2DXMatS2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      Mat2DXMatS2D%YX = M1%YX * M2%XX + M1%YY * M2%XY
      Mat2DXMatS2D%YY = M1%YX * M2%XY + M1%YY * M2%YY
   end function Mat2DXMatS2D

   function MatS2DXMatS2D(M1, M2)
      type(MatS2D), intent(IN)                     :: M1, M2
      type(Mat2D)                                  :: MatS2DXMatS2D

      MatS2DXMatS2D%XX = M1%XX * M2%XX + M1%XY * M2%XY
      MatS2DXMatS2D%XY = M1%XX * M2%XY + M1%XY * M2%YY
      MatS2DXMatS2D%YX = M1%XY * M2%XX + M1%YY * M2%XY
      MatS2DXMatS2D%YY = M1%XY * M2%XY + M1%YY * M2%YY
   end function MatS2DXMatS2D

   function Mat3DXMat3D(M1, M2)
      type(Mat3D), intent(IN)                      :: M1, M2
      type(Mat3D)                                  :: Mat3DXMat3D

      Mat3DXMat3D%XX = M1%XX * M2%XX + M1%XY * M2%YX + M1%XZ * M2%ZX
      Mat3DXMat3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%ZY
      Mat3DXMat3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ
      Mat3DXMat3D%YX = M1%YX * M2%XX + M1%YY * M2%YX + M1%YZ * M2%ZX
      Mat3DXMat3D%YY = M1%YX * M2%XY + M1%YY * M2%YY + M1%YZ * M2%ZY
      Mat3DXMat3D%YZ = M1%YX * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      Mat3DXMat3D%ZX = M1%ZX * M2%XX + M1%ZY * M2%YX + M1%ZZ * M2%ZX
      Mat3DXMat3D%ZY = M1%ZX * M2%XY + M1%ZY * M2%YY + M1%ZZ * M2%ZY
      Mat3DXMat3D%ZZ = M1%ZX * M2%XZ + M1%ZY * M2%YZ + M1%ZZ * M2%ZZ
   end function Mat3DXMat3D

   function Mat3DXMatS3D(M1, M2)
      type(Mat3D), intent(IN)                      :: M1
      type(MatS3D), intent(IN)                     :: M2
      type(Mat3D)                                  :: Mat3DXMatS3D

      Mat3DXMatS3D%XX = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ
      Mat3DXMatS3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%YZ
      Mat3DXMatS3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ
      Mat3DXMatS3D%YX = M1%YX * M2%XX + M1%YY * M2%XY + M1%YZ * M2%XZ
      Mat3DXMatS3D%YY = M1%YX * M2%XY + M1%YY * M2%YY + M1%YZ * M2%YZ
      Mat3DXMatS3D%YZ = M1%YX * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      Mat3DXMatS3D%ZX = M1%ZX * M2%XX + M1%ZY * M2%XY + M1%ZZ * M2%XZ
      Mat3DXMatS3D%ZY = M1%ZX * M2%XY + M1%ZY * M2%YY + M1%ZZ * M2%YZ
      Mat3DXMatS3D%ZZ = M1%ZX * M2%XZ + M1%ZY * M2%YZ + M1%ZZ * M2%ZZ
   end function Mat3DXMatS3D

   function MatS3DXMatS3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1, M2
      type(Mat3D)                                  :: MatS3DXMatS3D

      MatS3DXMatS3D%XX = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ
      MatS3DXMatS3D%XY = M1%XX * M2%XY + M1%XY * M2%YY + M1%XZ * M2%YZ
      MatS3DXMatS3D%XZ = M1%XX * M2%XZ + M1%XY * M2%YZ + M1%XZ * M2%ZZ
      MatS3DXMatS3D%YX = M1%XY * M2%XX + M1%YY * M2%XY + M1%YZ * M2%XZ
      MatS3DXMatS3D%YY = M1%XY * M2%XY + M1%YY * M2%YY + M1%YZ * M2%YZ
      MatS3DXMatS3D%YZ = M1%XY * M2%XZ + M1%YY * M2%YZ + M1%YZ * M2%ZZ
      MatS3DXMatS3D%ZX = M1%XZ * M2%XX + M1%YZ * M2%XY + M1%ZZ * M2%XZ
      MatS3DXMatS3D%ZY = M1%XZ * M2%XY + M1%YZ * M2%YY + M1%ZZ * M2%YZ
      MatS3DXMatS3D%ZZ = M1%XZ * M2%XZ + M1%YZ * M2%YZ + M1%ZZ * M2%ZZ
   end function MatS3DXMatS3D

   ! Overloading "/"
   function Vect2DQuot(V1, D1)
      type(Vect2D), intent(IN)                     :: V1
      PetscReal, intent(IN)                        :: D1
      type(Vect2D)                                 :: Vect2DQuot

      Vect2DQuot%X = V1%X / D1
      Vect2DQuot%Y = V1%Y / D1
   end function Vect2DQuot

   function Vect3DQuot(V1, D1)
      type(Vect3D), intent(IN)                     :: V1
      PetscReal, intent(IN)                        :: D1
      type(Vect3D)                                 :: Vect3DQuot

      Vect3DQuot%X = V1%X / D1
      Vect3DQuot%Y = V1%Y / D1
      Vect3DQuot%Z = V1%Z / D1
   end function Vect3DQuot

   function Mat2DQuot(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Mat2D), intent(IN)                      :: M1
      type(Mat2D)                                  :: Mat2DQuot

      Mat2DQuot%XX = M1%XX / D1
      Mat2DQuot%XY = M1%XY / D1
      Mat2DQuot%YX = M1%YX / D1
      Mat2DQuot%YY = M1%YY / D1
   end function Mat2DQuot

   function MatS2DQuot(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: MatS2DQuot

      MatS2DQuot%XX = M1%XX / D1
      MatS2DQuot%YY = M1%YY / D1
      MatS2DQuot%XY = M1%XY / D1
   end function MatS2DQuot

   function Mat3DQuot(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: Mat3DQuot

      Mat3DQuot%XX = M1%XX / D1
      Mat3DQuot%XY = M1%XY / D1
      Mat3DQuot%XZ = M1%XZ / D1
      Mat3DQuot%YX = M1%YX / D1
      Mat3DQuot%YY = M1%YY / D1
      Mat3DQuot%YZ = M1%YZ / D1
      Mat3DQuot%ZX = M1%ZX / D1
      Mat3DQuot%ZY = M1%ZY / D1
      Mat3DQuot%ZZ = M1%ZZ / D1
   end function Mat3DQuot

   function MatS3DQuot(M1, D1)
      PetscReal, intent(IN)                        :: D1
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: MatS3DQuot

      MatS3DQuot%XX = M1%XX / D1
      MatS3DQuot%YY = M1%YY / D1
      MatS3DQuot%ZZ = M1%ZZ / D1
      MatS3DQuot%YZ = M1%YZ / D1
      MatS3DQuot%XZ = M1%XZ / D1
      MatS3DQuot%XY = M1%XY / D1
   end function MatS3DQuot

   function Tens4OS2DQuot(T1, D1)
      type(Tens4OS2D), intent(IN)                  :: T1
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS2D)                              :: Tens4OS2DQuot

      Tens4OS2DQuot%XXXX = T1%XXXX / D1
      Tens4OS2DQuot%XXXY = T1%XXXY / D1
      Tens4OS2DQuot%XXYY = T1%XXYY / D1

      Tens4OS2DQuot%XYXY = T1%XYXY / D1
      Tens4OS2DQuot%YYXY = T1%YYXY / D1

      Tens4OS2DQuot%YYYY = T1%YYYY / D1
   end function Tens4OS2DQuot

   function Tens4OS3DQuot(T1, D1)
      type(Tens4OS3D), intent(IN)                  :: T1
      PetscReal, intent(IN)                        :: D1
      type(Tens4OS3D)                              :: Tens4OS3DQuot

      Tens4OS3DQuot%XXXX = T1%XXXX / D1
      Tens4OS3DQuot%XXYY = T1%XXYY / D1
      Tens4OS3DQuot%XXZZ = T1%XXZZ / D1
      Tens4OS3DQuot%XXYZ = T1%XXYZ / D1
      Tens4OS3DQuot%XXXZ = T1%XXXZ / D1
      Tens4OS3DQuot%XXXY = T1%XXXY / D1

      Tens4OS3DQuot%YYYY = T1%YYYY / D1
      Tens4OS3DQuot%YYZZ = T1%YYZZ / D1
      Tens4OS3DQuot%YYYZ = T1%YYYZ / D1
      Tens4OS3DQuot%YYXZ = T1%YYXZ / D1
      Tens4OS3DQuot%YYXY = T1%YYXY / D1

      Tens4OS3DQuot%ZZZZ = T1%ZZZZ / D1
      Tens4OS3DQuot%ZZYZ = T1%ZZYZ / D1
      Tens4OS3DQuot%ZZXZ = T1%ZZXZ / D1
      Tens4OS3DQuot%ZZXY = T1%ZZXY / D1

      Tens4OS3DQuot%YZYZ = T1%YZYZ / D1
      Tens4OS3DQuot%YZXZ = T1%YZXZ / D1
      Tens4OS3DQuot%YZXY = T1%YZXY / D1

      Tens4OS3DQuot%XZXZ = T1%XZXZ / D1
      Tens4OS3DQuot%XZXY = T1%XZXY / D1

      Tens4OS3DQuot%XYXY = T1%XYXY / D1
   end function Tens4OS3DQuot

   ! dot product in 2D and 3D
   function DotP2D(V1, V2)
      type(Vect2D), intent(IN)                     :: V1, V2
      PetscReal                                    :: DotP2D

      DotP2D = V1%X * V2%X + V1%Y * V2%Y
   end function DotP2D

   function DotP3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      PetscReal                                    :: DotP3D

      DotP3D = V1%X * V2%X + V1%Y * V2%Y + V1%Z * V2%Z
   end function DotP3D

   function ContP2D(M1, M2)
      ! tr(A^t x B)
      type(Mat2D), intent(IN)                      :: M1, M2
      PetscReal                                    :: ContP2D

      ContP2D = M1%XX * M2%XX + M1%XY * M2%XY + M1%YX * M2%YX + M1%YY * M2%YY
   end function ContP2D

   function ContP2DS(M1, M2)
      ! tr(A^t x B)
      type(MatS2D), intent(IN)                     :: M1, M2
      PetscReal                                    :: ContP2DS

      ContP2DS = M1%XX * M2%XX + M1%YY * M2%YY + 2.0_kr * M1%XY * M2%XY
   end function ContP2DS

#undef __FUNCT__
#define __FUNCT__ "Mat2DDotMatS2D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  Mat2DDotMatS2D:
!!!
   function Mat2DDotMatS2D(M1, M2)
      type(Mat2D), intent(IN)                      :: M1
      type(MatS2D), intent(IN)                     :: M2
      PetscReal                                    :: Mat2DDotMatS2D

      Mat2DDotMatS2D = symmetrize(M1) .DotP.M2
   end function Mat2DDotMatS2D

#undef __FUNCT__
#define __FUNCT__ "MatS2DDotMat2D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  MatS2DDotMat2D:
!!!
   function MatS2DDotMat2D(M1, M2)
      type(MatS2D), intent(IN)                     :: M1
      type(Mat2D), intent(IN)                      :: M2
      PetscReal                                    :: MatS2DDotMat2D

      MatS2DDotMat2D = M1.DotP.symmetrize(M2)
   end function MatS2DDotMat2D

#undef __FUNCT__
#define __FUNCT__ "Mat2DDotMatS3D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  Mat3DDotMatS3D:
!!!
   function Mat3DDotMatS3D(M1, M2)
      type(Mat3D), intent(IN)                      :: M1
      type(MatS3D), intent(IN)                     :: M2
      PetscReal                                    :: Mat3DDotMatS3D

      Mat3DDotMatS3D = symmetrize(M1) .DotP.M2
   end function Mat3DDotMatS3D

#undef __FUNCT__
#define __FUNCT__ "MatS3DDotMat3D"
!!! author: Blaise Bourdin (2016, bourdin@lsu.edu)
!!!
!!!  MatS3DDotMat3D:
!!!
   function MatS3DDotMat3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1
      type(Mat3D), intent(IN)                      :: M2
      PetscReal                                    :: MatS3DDotMat3D

      MatS3DDotMat3D = M1.DotP.symmetrize(M2)
   end function MatS3DDotMat3D

   function ContP3D(M1, M2)
      ! tr(A^t x B)
      type(Mat3D), intent(IN)                      :: M1, M2
      PetscReal                                   :: ContP3D

      ContP3D = M1%XX * M2%XX + M1%XY * M2%XY + M1%XZ * M2%XZ + &
           &    M1%YX * M2%YX + M1%YY * M2%YY + M1%YZ * M2%YZ + &
           &    M1%ZX * M2%ZX + M1%ZY * M2%ZY + M1%ZZ * M2%ZZ
   end function ContP3D

   function ContP3DS(M1, M2)
      ! tr(A^t x B)
      type(MatS3D), intent(IN)                     :: M1, M2
      PetscReal                                    :: ContP3DS

      ContP3DS = M1%XX * M2%XX + M1%YY * M2%YY + M1%ZZ * M2%ZZ + 2.0_kr * M1%YZ * M2%YZ &
                 + 2.0_kr * M1%XZ * M2%XZ + 2.0_kr * M1%XY * M2%XY
   end function ContP3DS

   ! cross product 3D
   function CrossP3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      type(Vect3D)                                 :: CrossP3D

      CrossP3D%X = V1%Y * V2%Z - V1%Z * V2%Y
      CrossP3D%Y = V1%Z * V2%X - V1%X * V2%Z
      CrossP3D%Z = V1%X * V2%Y - V1%Y * V2%X
   end function CrossP3D

   ! Transpose
   function Transpose2D(M1)
      type(Mat2D), intent(IN)                     :: M1
      type(Mat2D)                                 :: Transpose2D

      Transpose2D%XX = M1%XX
      Transpose2D%XY = M1%YX
      Transpose2D%YX = M1%XY
      Transpose2D%YY = M1%YY
   end function Transpose2D

   function Transpose3D(M1)
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: Transpose3D

      Transpose3D%XX = M1%XX
      Transpose3D%XY = M1%YX
      Transpose3D%XZ = M1%ZX
      Transpose3D%YX = M1%XY
      Transpose3D%YY = M1%YY
      Transpose3D%YZ = M1%ZY
      Transpose3D%ZX = M1%XZ
      Transpose3D%ZY = M1%YZ
      Transpose3D%ZZ = M1%ZZ
   end function Transpose3D

   ! Tensor product
   function TensPVect2D(V1, V2)
      type(Vect2D), intent(IN)                     :: V1
      type(Vect2D), intent(IN)                     :: V2
      type(Mat2D)                                  :: TensPVect2D

      TensPVect2D%XX = V1%X * V2%X
      TensPVect2D%XY = V1%X * V2%Y
      TensPVect2D%YX = V1%Y * V2%X
      TensPVect2D%YY = V1%Y * V2%Y
   end function TensPVect2D

   function TensPVect3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      type(Mat3D)                                  :: TensPVect3D

      TensPVect3D%XX = V1%X * V2%X
      TensPVect3D%XY = V1%X * V2%Y
      TensPVect3D%XZ = V1%X * V2%Z
      TensPVect3D%YX = V1%Y * V2%X
      TensPVect3D%YY = V1%Y * V2%Y
      TensPVect3D%YZ = V1%Y * V2%Z
      TensPVect3D%ZX = V1%Z * V2%X
      TensPVect3D%ZY = V1%Z * V2%Y
      TensPVect3D%ZZ = V1%Z * V2%Z
   end function TensPVect3D

   ! Symmetrized product
   function SymPVect2D(V1, V2)
      type(Vect2D), intent(IN)                     :: V1, V2
      type(MatS2D)                                 :: SymPVect2D

      SymPVect2D = Symmetrize(V1.TensP.V2)
   end function SymPVect2D

   function SymPVect3D(V1, V2)
      type(Vect3D), intent(IN)                     :: V1, V2
      type(MatS3D)                                 :: SymPVect3D

      SymPVect3D = Symmetrize(V1.TensP.V2)
   end function SymPVect3D

   function SymPMatS2D(M1, M2)
      type(MatS2D), intent(IN)                     :: M1, M2
      type(Tens4OS2D)                              :: SymPMatS2D

      SymPMatS2D%XXXX = M1%XX * M2%XX
      SymPMatS2D%XXYY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_kr
      SymPMatS2D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_kr

      SymPMatS2D%YYYY = M1%YY * M2%YY
      SymPMatS2D%YYXY = (M1%YY * M2%XY + M1%XY * M2%YY) * 0.5_kr

      SymPMatS2D%XYXY = M1%XY * M2%XY
   end function SymPMatS2D

   function SymPMatS3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1, M2
      type(Tens4OS3D)                              :: SymPMatS3D

      SymPMatS3D%XXXX = M1%XX * M2%XX
      SymPMatS3D%XXYY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_kr
      SymPMatS3D%XXZZ = (M1%XX * M2%ZZ + M1%ZZ * M2%XX) * 0.5_kr
      SymPMatS3D%XXYZ = (M1%XX * M2%YZ + M1%YZ * M2%XX) * 0.5_kr
      SymPMatS3D%XXXZ = (M1%XX * M2%XZ + M1%XZ * M2%XX) * 0.5_kr
      SymPMatS3D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_kr

      SymPMatS3D%YYYY = M1%YY * M2%YY
      SymPMatS3D%YYZZ = (M1%YY * M2%ZZ + M1%ZZ * M2%YY) * 0.5_kr
      SymPMatS3D%YYYZ = (M1%YY * M2%YZ + M1%YZ * M2%YY) * 0.5_kr
      SymPMatS3D%YYXZ = (M1%YY * M2%XZ + M1%XZ * M2%YY) * 0.5_kr
      SymPMatS3D%YYXY = (M1%YY * M2%XY + M1%XY * M2%YY) * 0.5_kr

      SymPMatS3D%ZZZZ = M1%ZZ * M2%ZZ
      SymPMatS3D%ZZYZ = (M1%ZZ * M2%YZ + M1%YZ * M2%ZZ) * 0.5_kr
      SymPMatS3D%ZZXZ = (M1%ZZ * M2%XZ + M1%XZ * M2%ZZ) * 0.5_kr
      SymPMatS3D%ZZXY = (M1%ZZ * M2%XY + M1%XY * M2%ZZ) * 0.5_kr

      SymPMatS3D%YZYZ = M1%YZ * M2%YZ
      SymPMatS3D%YZXZ = (M1%YZ * M2%XZ + M1%XZ * M2%YZ) * 0.5_kr
      SymPMatS3D%YZXY = (M1%YZ * M2%XY + M1%XY * M2%YZ) * 0.5_kr

      SymPMatS3D%XZXZ = M1%XZ * M2%XZ
      SymPMatS3D%XZXY = (M1%XZ * M2%XY + M1%XY * M2%XZ) * 0.5_kr

      SymPMatS3D%XYXY = M1%XY * M2%XY
   end function SymPMatS3D

   function oDotMatS2D(M1, M2)
      type(MatS2D), intent(IN)                     :: M1, M2
      type(Tens4OS2D)                              :: oDotMatS2D

      oDotMatS2D%XXXX = M1%XX * M2%XX
      oDotMatS2D%XXYY = M1%XY * M2%XY
      oDotMatS2D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_kr

      oDotMatS2D%YYYY = M1%YY * M2%YY
      oDotMatS2D%YYXY = (M1%XY * M2%YY + M1%YY * M2%XY) * 0.5_kr

      oDotMatS2D%XYXY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_kr
   end function oDotMatS2D

   function oDotMatS3D(M1, M2)
      type(MatS3D), intent(IN)                     :: M1, M2
      type(Tens4OS3D)                              :: oDotMatS3D

      oDotMatS3D%XXXX = M1%XX * M2%XX
      oDotMatS3D%XXYY = M1%XY * M2%XY
      oDotMatS3D%XXZZ = M1%XZ * M2%XZ
      oDotMatS3D%XXYZ = (M1%XY * M2%XZ + M1%XZ * M2%YZ) * 0.5_kr
      oDotMatS3D%XXXZ = (M1%XX * M2%XZ + M1%XZ * M2%XX) * 0.5_kr
      oDotMatS3D%XXXY = (M1%XX * M2%XY + M1%XY * M2%XX) * 0.5_kr

      oDotMatS3D%YYYY = M1%YY * M2%YY
      oDotMatS3D%YYZZ = M1%YZ * M2%YZ
      oDotMatS3D%YYYZ = (M1%YY * M2%YZ + M1%YZ * M2%YY) * 0.5_kr
      oDotMatS3D%YYXZ = (M1%XY * M2%YZ + M1%YZ * M2%XY) * 0.5_kr
      oDotMatS3D%YYXY = (M1%XY * M2%YY + M1%YY * M2%XY) * 0.5_kr

      oDotMatS3D%ZZZZ = M1%ZZ * M2%ZZ
      oDotMatS3D%ZZYZ = (M1%YZ * M2%ZZ + M1%ZZ * M2%YZ) * 0.5_kr
      oDotMatS3D%ZZXZ = (M1%XZ * M2%ZZ + M1%ZZ * M2%XZ) * 0.5_kr
      oDotMatS3D%ZZXY = (M1%XZ * M2%YZ + M1%YZ * M2%XZ) * 0.5_kr

      oDotMatS3D%YZYZ = (M1%YY * M2%ZZ + M1%ZZ * M2%YY) * 0.5_kr
      oDotMatS3D%YZXZ = (M1%XY * M2%ZZ + M1%ZZ * M2%XY) * 0.5_kr
      oDotMatS3D%YZXY = (M1%XY * M2%YZ + M1%YZ * M2%XY) * 0.5_kr

      oDotMatS3D%XZXZ = (M1%XX * M2%ZZ + M1%ZZ * M2%XX) * 0.5_kr
      oDotMatS3D%XZXY = (M1%XX * M2%YZ + M1%YZ * M2%XX) * 0.5_kr

      oDotMatS3D%XYXY = (M1%XX * M2%YY + M1%YY * M2%XX) * 0.5_kr
   end function oDotMatS3D

   function Trace2D(M1)
      type(Mat2D), intent(IN)                      :: M1
      PetscReal                                    :: Trace2D

      Trace2D = M1%XX + M1%YY
   end function Trace2D

   function Trace2DS(M1)
      type(MatS2D), intent(IN)                     :: M1
      PetscReal                                    :: Trace2DS

      Trace2DS = M1%XX + M1%YY
   end function Trace2DS

   function Trace3D(M1)
      type(Mat3D), intent(IN)                      :: M1
      PetscReal                                    :: Trace3D

      Trace3D = M1%XX + M1%YY + M1%ZZ
   end function Trace3D

   function Trace3DS(M1)
      type(MatS3D), intent(IN)                     :: M1
      PetscReal                                    :: Trace3DS

      Trace3DS = M1%XX + M1%YY + M1%ZZ
   end function Trace3DS

   subroutine Vect2D_Get_Real(V1, R1)
      type(Vect2D), intent(OUT)                    :: V1
      PetscReal, intent(IN)                        :: R1

      V1%X = R1
      V1%Y = R1
   end subroutine Vect2D_Get_Real

   subroutine Vect3D_Get_Real(V1, R1)
      type(Vect3D), intent(OUT)                    :: V1
      PetscReal, intent(IN)                        :: R1

      V1%X = R1
      V1%Y = R1
      V1%Z = R1
   end subroutine Vect3D_Get_Real

   subroutine Vect2D_Get_VectR(V1, R1)
      type(Vect2D), intent(OUT)               :: V1
      PetscReal, dimension(2), intent(IN)     :: R1

      V1%X = R1(1)
      V1%Y = R1(2)
   end subroutine Vect2D_Get_VectR

   subroutine Vect3D_Get_VectR(V1, R1)
      type(Vect3D), intent(OUT)               :: V1
      PetscReal, dimension(3), intent(IN)     :: R1

      V1%X = R1(1)
      V1%Y = R1(2)
      V1%Z = R1(3)
   end subroutine Vect3D_Get_VectR

   subroutine Vect2DEQ(V1, V2)
      type(Vect2D), intent(OUT)                    :: V1
      type(Vect2D), intent(IN)                     :: V2

      V1%X = V2%X
      V1%Y = V2%Y
   end subroutine Vect2DEQ

   subroutine Vect3DEQ(V1, V2)
      type(Vect3D), intent(OUT)                    :: V1
      type(Vect3D), intent(IN)                     :: V2

      V1%X = V2%X
      V1%Y = V2%Y
      V1%Z = V2%Z
   end subroutine Vect3DEQ

   subroutine Mat2D_Get_Real(M1, R1)
      type(Mat2D), intent(OUT)                     :: M1
      PetscReal, intent(IN)                        :: R1

      M1%XX = R1; M1%XY = R1
      M1%YX = R1; M1%YY = R1
   end subroutine Mat2D_Get_Real

   subroutine Mat3D_Get_Real(M1, R1)
      type(Mat3D), intent(OUT)                     :: M1
      PetscReal, intent(IN)                        :: R1

      M1%XX = R1; M1%XY = R1; M1%XZ = R1
      M1%YX = R1; M1%YY = R1; M1%YZ = R1
      M1%ZX = R1; M1%ZY = R1; M1%ZZ = R1
   end subroutine Mat3D_Get_Real

   subroutine Mat2DEQ(M1, M2)
      type(Mat2D), intent(OUT)                     :: M1
      type(Mat2D), intent(IN)                      :: M2

      M1%XX = M2%XX; M1%XY = M2%XY
      M1%YX = M2%YX; M1%YY = M2%YY
   end subroutine Mat2DEQ

   subroutine Mat3DEQ(M1, M2)
      type(Mat3D), intent(OUT)                     :: M1
      type(Mat3D), intent(IN)                      :: M2

      M1%XX = M2%XX; M1%XY = M2%XY; M1%XZ = M2%XZ
      M1%YX = M2%YX; M1%YY = M2%YY; M1%YZ = M2%YZ
      M1%ZX = M2%ZX; M1%ZY = M2%ZY; M1%ZZ = M2%ZZ
   end subroutine Mat3DEQ

   subroutine MatS2D_Get_Real(M1, R1)
      type(MatS2D), intent(OUT)                    :: M1
      PetscReal, intent(IN)                        :: R1

      M1%XX = R1; M1%YY = R1; M1%XY = R1
   end subroutine MatS2D_Get_Real

   subroutine MatS3D_Get_Real(M1, R1)
      type(MatS3D), intent(OUT)                    :: M1
      PetscReal, intent(IN)                        :: R1

      M1%XX = R1; M1%YY = R1; M1%ZZ = R1
      M1%YZ = R1; M1%XZ = R1; M1%XY = R1
   end subroutine MatS3D_Get_Real

   subroutine MatS2DEQ(M1, M2)
      type(MatS2D), intent(OUT)                    :: M1
      type(MatS2D), intent(IN)                     :: M2

      M1%XX = M2%XX; M1%YY = M2%YY; M1%XY = M2%XY
   end subroutine MatS2DEQ

   subroutine MatS3DEQ(M1, M2)
      type(MatS3D), intent(OUT)                    :: M1
      type(MatS3D), intent(IN)                     :: M2

      M1%XX = M2%XX; M1%YY = M2%YY; M1%ZZ = M2%ZZ
      M1%YZ = M2%YZ; M1%XZ = M2%XZ; M1%XY = M2%XY
   end subroutine MatS3DEQ

   subroutine MatS2D_Get_VectR(M1, R1)
      type(MatS2D), intent(OUT)                    :: M1
      PetscReal, dimension(3), intent(IN)          :: R1

      M1%XX = R1(1)
      M1%YY = R1(2)
      M1%XY = R1(3)
   end subroutine MatS2D_Get_VectR

   subroutine MatS3D_Get_VectR(M1, R1)
      type(MatS3D), intent(OUT)                    :: M1
      PetscReal, dimension(6), intent(IN)          :: R1

      M1%XX = R1(1)
      M1%YY = R1(2)
      M1%ZZ = R1(3)
      M1%YZ = R1(4)
      M1%XZ = R1(5)
      M1%XY = R1(6)
   end subroutine MatS3D_Get_VectR

   subroutine VectR_Get_MatS2D(R1, M1)
      PetscReal, dimension(3), intent(OUT)          :: R1
      type(MatS2D), intent(IN)                      :: M1

      R1(1) = M1%XX
      R1(2) = M1%YY
      R1(3) = M1%XY
   end subroutine VectR_Get_MatS2D

   subroutine VectR_Get_MatS3D(R1, M1)
      PetscReal, dimension(6), intent(OUT)          :: R1
      type(MatS3D), intent(IN)                      :: M1

      R1(1) = M1%XX
      R1(2) = M1%YY
      R1(3) = M1%ZZ
      R1(4) = M1%YZ
      R1(5) = M1%XZ
      R1(6) = M1%XY
   end subroutine VectR_Get_MatS3D

   subroutine Mat2DGetArray(M, A)
      type(Mat2D), intent(OUT)                     :: M
      PetscReal, dimension(2, 2), intent(IN)       :: A

      M%XX = A(1, 1)
      M%XY = A(1, 2)
      M%YX = A(2, 1)
      M%YY = A(2, 2)
   end subroutine Mat2DGetArray

   subroutine MatS2DGetArray(M, A)
      type(MatS2D), intent(OUT)                    :: M
      PetscReal, dimension(2, 2), intent(IN)       :: A

      M%XX = A(1, 1)
      M%YY = A(2, 2)
      M%XY = A(1, 2)
   end subroutine MatS2DGetArray

   subroutine Mat3DGetArray(M, A)
      type(Mat3D), intent(OUT)                     :: M
      PetscReal, dimension(3, 3), intent(IN)       :: A

      M%XX = A(1, 1)
      M%XY = A(1, 2)
      M%XZ = A(1, 3)
      M%YX = A(2, 1)
      M%YY = A(2, 2)
      M%YZ = A(2, 3)
      M%ZX = A(3, 1)
      M%ZY = A(3, 2)
      M%ZZ = A(3, 3)
   end subroutine Mat3DGetArray

   subroutine MatS3DGetArray(M, A)
      type(MatS3D), intent(OUT)                    :: M
      PetscReal, dimension(3, 3), intent(IN)       :: A

      M%XX = A(1, 1)
      M%YY = A(2, 2)
      M%ZZ = A(3, 3)
      M%YZ = A(2, 3)
      M%XZ = A(1, 3)
      M%XY = A(1, 2)
   end subroutine MatS3DGetArray

   subroutine ArrayGetMat2D(A, M)
      PetscReal, dimension(2, 2), intent(OUT)        :: A
      type(Mat2D), intent(IN)                        :: M

      A(1, 1) = M%XX
      A(1, 2) = M%XY
      A(2, 1) = M%YX
      A(2, 2) = M%YY
   end subroutine ArrayGetMat2D

   subroutine ArrayGetMatS2D(A, M)
      PetscReal, dimension(2, 2), intent(OUT)        :: A
      type(MatS2D), intent(IN)                       :: M

      A(1, 1) = M%XX
      A(1, 2) = M%XY
      A(2, 1) = M%XY
      A(2, 2) = M%YY
   end subroutine ArrayGetMatS2D

   subroutine ArrayGetMat3D(A, M)
      PetscReal, dimension(3, 3), intent(OUT)        :: A
      type(Mat3D), intent(IN)                        :: M

      A(1, 1) = M%XX
      A(1, 2) = M%XY
      A(1, 3) = M%XZ
      A(2, 1) = M%YX
      A(2, 2) = M%YY
      A(2, 3) = M%YZ
      A(3, 1) = M%ZX
      A(3, 2) = M%ZY
      A(3, 3) = M%ZZ
   end subroutine ArrayGetMat3D

   subroutine ArrayGetMatS3D(A, M)
      PetscReal, dimension(3, 3), intent(OUT)        :: A
      type(MatS3D), intent(IN)                       :: M

      A(1, 1) = M%XX
      A(1, 2) = M%XY
      A(1, 3) = M%XZ
      A(2, 1) = M%XY
      A(2, 2) = M%YY
      A(2, 3) = M%YZ
      A(3, 1) = M%XZ
      A(3, 2) = M%YZ
      A(3, 3) = M%ZZ
   end subroutine ArrayGetMatS3D

   subroutine Tens4OS2DEQ(T1, T2)
      type(Tens4OS2D), intent(OUT)                 :: T1
      type(Tens4OS2D), intent(IN)                  :: T2

      T1%XXXX = T2%XXXX
      T1%XXXY = T2%XXXY
      T1%XXYY = T2%XXYY

      T1%XYXY = T2%XYXY
      T1%YYXY = T2%YYXY

      T1%YYYY = T2%YYYY
   end subroutine Tens4OS2DEQ

   subroutine Tens4OS3DEQ(T1, T2)
      type(Tens4OS3D), intent(OUT)                 :: T1
      type(Tens4OS3D), intent(IN)                  :: T2

      T1%XXXX = T2%XXXX
      T1%XXYY = T2%XXYY
      T1%XXZZ = T2%XXZZ
      T1%XXYZ = T2%XXYZ
      T1%XXXZ = T2%XXXZ
      T1%XXXY = T2%XXXY

      T1%YYYY = T2%YYYY
      T1%YYZZ = T2%YYZZ
      T1%YYYZ = T2%YYYZ
      T1%YYXZ = T2%YYXZ
      T1%YYXY = T2%YYXY

      T1%ZZZZ = T2%ZZZZ
      T1%ZZYZ = T2%ZZYZ
      T1%ZZXZ = T2%ZZXZ
      T1%ZZXY = T2%ZZXY

      T1%YZYZ = T2%YZYZ
      T1%YZXZ = T2%YZXZ
      T1%YZXY = T2%YZXY

      T1%XZXZ = T2%XZXZ
      T1%XZXY = T2%XZXY

      T1%XYXY = T2%XYXY
   end subroutine Tens4OS3DEQ

   subroutine Tens4OS2D_Get_Real(T1, D1)
      type(Tens4OS2D), intent(OUT)                 :: T1
      PetscReal, intent(IN)                        :: D1

      T1%XXXX = D1
      T1%XXXY = D1
      T1%XXYY = D1

      T1%XYXY = D1
      T1%YYXY = D1

      T1%YYYY = D1
   end subroutine Tens4OS2D_Get_Real

   subroutine Tens4OS3D_Get_Real(T1, D1)
      type(Tens4OS3D), intent(OUT)                 :: T1
      PetscReal, intent(IN)                        :: D1

      T1%XXXX = D1
      T1%XXYY = D1
      T1%XXZZ = D1
      T1%XXYZ = D1
      T1%XXXZ = D1
      T1%XXXY = D1

      T1%YYYY = D1
      T1%YYZZ = D1
      T1%YYYZ = D1
      T1%YYXZ = D1
      T1%YYXY = D1

      T1%ZZZZ = D1
      T1%ZZYZ = D1
      T1%ZZXZ = D1
      T1%ZZXY = D1

      T1%YZYZ = D1
      T1%YZXZ = D1
      T1%YZXY = D1

      T1%XZXZ = D1
      T1%XZXY = D1

      T1%XYXY = D1
   end subroutine Tens4OS3D_Get_Real

   subroutine Tens4OS2DToArray(A, T)
      PetscReal, dimension(3, 3), intent(OUT)        :: A
      type(Tens4OS2D), intent(IN)                    :: T

      A(1, 1) = T%XXXX
      A(1, 2) = T%XXYY
      A(1, 3) = T%XXXY * 2.0_kr

      A(2, 1) = A(1, 2)
      A(2, 2) = T%YYYY
      A(2, 3) = T%YYXY * 2.0_kr

      A(3, 1) = A(1, 3)
      A(3, 2) = A(2, 3)
      A(3, 3) = T%XYXY * 2.0_kr
   end subroutine Tens4OS2DToArray

   subroutine ArrayToTens4OS2D(T, A)
      type(Tens4OS2D), intent(OUT)                   :: T
      PetscReal, dimension(3, 3), intent(IN)         :: A

      T%XXXX = A(1, 1)
      T%XXYY = A(1, 2)
      T%XXXY = A(1, 3)*.5_kr

      T%YYYY = A(2, 2)
      T%YYXY = A(2, 3)*.5_kr

      T%XYXY = A(3, 3)*.5_kr
   end subroutine ArrayToTens4OS2D

   subroutine Tens4OS3DToArray(A, T)
      PetscReal, dimension(6, 6), intent(OUT)        :: A
      type(Tens4OS3D), intent(IN)                    :: T

      A(1, 1) = T%XXXX
      A(1, 2) = T%XXYY
      A(1, 3) = T%XXZZ
      A(1, 4) = T%XXYZ * 2.0_kr
      A(1, 5) = T%XXXZ * 2.0_kr
      A(1, 6) = T%XXXY * 2.0_kr

      A(2, 1) = A(1, 2)
      A(2, 2) = T%YYYY
      A(2, 3) = T%YYZZ
      A(2, 4) = T%YYYZ * 2.0_kr
      A(2, 5) = T%YYXZ * 2.0_kr
      A(2, 6) = T%YYXY * 2.0_kr

      A(3, 1) = A(1, 3)
      A(3, 2) = A(2, 3)
      A(3, 3) = T%ZZZZ
      A(3, 4) = T%ZZYZ * 2.0_kr
      A(3, 5) = T%ZZXZ * 2.0_kr
      A(3, 6) = T%ZZXY * 2.0_kr

      A(4, 1) = A(1, 4)
      A(4, 2) = A(2, 4)
      A(4, 3) = A(3, 4)
      A(4, 4) = T%YZYZ * 2.0_kr
      A(4, 5) = T%YZXZ * 2.0_kr
      A(4, 6) = T%YZXY * 2.0_kr

      A(5, 1) = A(1, 5)
      A(5, 2) = A(2, 5)
      A(5, 3) = A(3, 5)
      A(5, 4) = A(4, 5)
      A(5, 5) = T%XZXZ * 2.0_kr
      A(5, 6) = T%XZXY * 2.0_kr

      A(6, 1) = A(1, 6)
      A(6, 2) = A(2, 6)
      A(6, 3) = A(3, 6)
      A(6, 4) = A(4, 6)
      A(6, 5) = A(5, 6)
      A(6, 6) = T%XYXY * 2.0_kr
   end subroutine Tens4OS3DToArray

   subroutine ArrayToTens4OS3D(T, A)
      type(Tens4OS3D), intent(OUT)                   :: T
      PetscReal, dimension(6, 6), intent(IN)         :: A

      T%XXXX = A(1, 1)
      T%XXYY = A(1, 2)
      T%XXZZ = A(1, 3)
      T%XXYZ = A(1, 4)*.5_kr
      T%XXXZ = A(1, 5)*.5_kr
      T%XXXY = A(1, 6)*.5_kr

      T%YYYY = A(2, 2)
      T%YYZZ = A(2, 3)
      T%YYYZ = A(2, 4)*.5_kr
      T%YYXZ = A(2, 5)*.5_kr
      T%YYXY = A(2, 6)*.5_kr

      T%ZZZZ = A(3, 3)
      T%ZZYZ = A(3, 4)*.5_kr
      T%ZZXZ = A(3, 5)*.5_kr
      T%ZZXY = A(3, 6)*.5_kr

      T%YZYZ = A(4, 4)*.5_kr
      T%YZXZ = A(4, 5)*.5_kr
      T%YZXY = A(4, 6)*.5_kr

      T%XZXZ = A(5, 5)*.5_kr
      T%XZXY = A(5, 6)*.5_kr

      T%XYXY = A(6, 6)*.5_kr
   end subroutine ArrayToTens4OS3D

   !!! Overloading euclidian norm of derived types
   PetscReal function Vect2DNorm(V)
      type(Vect2D), intent(IN)                     :: V

      Vect2DNorm = sqrt(V%X**2 + V%Y**2)
   end function Vect2DNorm

   PetscReal function Vect3DNorm(V)
      type(Vect3D), intent(IN)                     :: V

      Vect3DNorm = sqrt(V%X**2 + V%Y**2 + V%Z**2)
   end function Vect3DNorm

   PetscReal function Mat2DNorm(M)
      type(Mat2D), intent(IN)                      :: M

      Mat2DNorm = sqrt(M%XX**2 + M%XY**2 + M%YX**2 + M%YY**2)
   end function Mat2DNorm

   PetscReal function MatS2DNorm(M)
      type(MatS2D), intent(IN)                     :: M

      MatS2DNorm = sqrt(M%XX**2 + 2.0_kr * M%XY**2 + M%YY**2)
   end function MatS2DNorm

   PetscReal function Mat3DNorm(M)
      type(Mat3D), intent(IN)                      :: M

      Mat3DNorm = sqrt(M%XX**2 + M%XY**2 + M%XZ**2 + M%YX**2 + M%YY**2 + M%YZ**2 + M%ZX**2 + M%ZY**2 + M%ZZ**2)
   end function Mat3DNorm

   PetscReal function MatS3DNorm(M)
      type(MatS3D), intent(IN)                     :: M

      MatS3DNorm = sqrt(M%XX**2 + 2.0_kr * M%XY**2 + 2.0_kr * M%XZ**2 + M%YY**2 + 2.0_kr * M%YZ**2 + M%ZZ**2)
   end function MatS3DNorm

   function Symmetrize2D(M1)
      type(Mat2D), intent(IN)                      :: M1
      type(MatS2D)                                 :: Symmetrize2D

      Symmetrize2D%XX = M1%XX
      Symmetrize2D%YY = M1%YY
      Symmetrize2D%XY = (M1%XY + M1%YX) * 0.5_kr
   end function Symmetrize2D

   function Symmetrize3D(M1)
      type(Mat3D), intent(IN)                      :: M1
      type(MatS3D)                                 :: Symmetrize3D

      Symmetrize3D%XX = M1%XX
      Symmetrize3D%YY = M1%YY
      Symmetrize3D%ZZ = M1%ZZ

      Symmetrize3D%YZ = (M1%YZ + M1%ZY) * 0.5_kr
      Symmetrize3D%XZ = (M1%XZ + M1%ZX) * 0.5_kr
      Symmetrize3D%XY = (M1%XY + M1%YX) * 0.5_kr
   end function Symmetrize3D

   subroutine MatS2DToMat2D(M1, M2)
      type(Mat2D), intent(OUT)                     :: M1
      type(MatS2D), intent(IN)                     :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%YX = M2%XY
      M1%YY = M2%YY
   end subroutine MatS2DToMat2D

   subroutine MatS3DToMat3D(M1, M2)
      type(Mat3D), intent(OUT)                     :: M1
      type(MatS3D), intent(IN)                     :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%XZ = M2%XZ
      M1%YX = M2%XY
      M1%YY = M2%YY
      M1%YZ = M2%YZ
      M1%ZX = M2%XZ
      M1%ZY = M2%YZ
      M1%ZZ = M2%ZZ
   end subroutine MatS3DToMat3D

   subroutine Mat2DToMatS2D(M1, M2)
      type(MatS2D), intent(OUT)                    :: M1
      type(Mat2D), intent(IN)                      :: M2

      M1%XX = M2%XX
      M1%XY = M2%XY
      M1%YY = M2%YY
   end subroutine Mat2DToMatS2D

   subroutine Mat3DToMatS3D(M1, M2)
      type(MatS3D), intent(OUT)                    :: M1
      type(Mat3D), intent(IN)                      :: M2

      M1%XX = M2%XX
      M1%YY = M2%YY
      M1%ZZ = M2%ZZ
      M1%YZ = M2%YZ
      M1%XZ = M2%XZ
      M1%XY = M2%XY
   end subroutine Mat3DToMatS3D

   function RARtMat2D(A, R)
      ! A <- R.A.R^T
      type(Mat2D), intent(IN)                      :: A
      type(Mat2D), intent(IN)                      :: R
      type(Mat2D)                                  :: RARtMat2D

      RARtMat2D = R * A * transpose(R)
   end function RaRtMat2D

   function RARtMatS2D(A, R)
      ! A <- R.A.R^T
      type(MatS2D), intent(IN)                     :: A
      type(Mat2D), intent(IN)                      :: R
      type(MatS2D)                                 :: RARtMatS2D

      RARtMatS2D = R * A * transpose(R)
   end function RaRtMatS2D

   function RARtMat3D(A, R)
      ! A <- R.A.R^T
      type(Mat3D), intent(IN)                      :: A
      type(Mat3D), intent(IN)                      :: R
      type(Mat3D)                                  :: RARtMat3D

      RARtMat3D = R * A * transpose(R)
   end function RaRtMat3D

   function RARtMatS3D(A, R)
      ! A <- R.A.R^T
      type(MatS3D), intent(IN)                     :: A
      type(Mat3D), intent(IN)                      :: R
      type(MatS3D)                                 :: RARtMatS3D

      RARtMatS3D = R * A * transpose(R)
   end function RaRtMatS3D

   function RtARMat2D(A, R)
      ! A <- Rt.A.R
      type(Mat2D), intent(IN)                      :: A
      type(Mat2D), intent(IN)                      :: R
      type(Mat2D)                                  :: RtARMat2D

      RtARMat2D = transpose(R) * A * R
   end function RtaRMat2D

   function RtaRMatS2D(A, R)
      ! A <- R^T.A.R
      type(MatS2D), intent(IN)                     :: A
      type(Mat2D), intent(IN)                      :: R
      type(MatS2D)                                 :: RtaRMatS2D

      RtaRMatS2D = transpose(R) * A * R
   end function RtaRMatS2D

   function RtaRMat3D(A, R)
      ! A <- R^T.A.R
      type(Mat3D), intent(IN)                      :: A
      type(Mat3D), intent(IN)                      :: R
      type(Mat3D)                                  :: RtaRMat3D

      RtaRMat3D = transpose(R) * A * R
   end function RtaRMat3D

   function RtaRMatS3D(A, R)
      ! A <- R^T.A.R
      type(MatS3D), intent(IN)                     :: A
      type(Mat3D), intent(IN)                      :: R
      type(MatS3D)                                 :: RtaRMatS3D

      RtaRMatS3D = transpose(R) * A * R
   end function RtaRMatS3D

   function DeviatoricPart2D(M1)
      type(Mat2D), intent(IN)                      :: M1
      type(Mat2D)                                  :: DeviatoricPart2D

      PetscReal                                   :: M1_Trace

      M1_Trace = Trace(M1)
      DeviatoricPart2D%XX = M1%XX - M1_Trace * 0.5_kr
      DeviatoricPart2D%XY = M1%XY
      DeviatoricPart2D%YX = M1%YX
      DeviatoricPart2D%YY = M1%YY - M1_Trace * 0.5_kr
   end function DeviatoricPart2D

   function DeviatoricPart2DS(M1)
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: DeviatoricPart2DS

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      DeviatoricPart2DS%XX = M1%XX - M1_Trace * 0.5_kr
      DeviatoricPart2DS%YY = M1%YY - M1_Trace * 0.5_kr
      DeviatoricPart2DS%XY = M1%XY
   end function DeviatoricPart2DS

   function DeviatoricPart3D(M1)
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: DeviatoricPart3D

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      DeviatoricPart3D%XX = M1%XX - M1_Trace / 3.0_kr
      DeviatoricPart3D%XY = M1%XY
      DeviatoricPart3D%XZ = M1%XZ

      DeviatoricPart3D%YX = M1%YX
      DeviatoricPart3D%YY = M1%YY - M1_Trace / 3.0_kr
      DeviatoricPart3D%YZ = M1%YZ

      DeviatoricPart3D%ZX = M1%ZX
      DeviatoricPart3D%ZY = M1%ZY
      DeviatoricPart3D%ZZ = M1%ZZ - M1_Trace / 3.0_kr
   end function DeviatoricPart3D

   function DeviatoricPart3DS(M1)
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: DeviatoricPart3DS

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      DeviatoricPart3DS%XX = M1%XX - M1_Trace / 3.0_kr
      DeviatoricPart3DS%YY = M1%YY - M1_Trace / 3.0_kr
      DeviatoricPart3DS%ZZ = M1%ZZ - M1_Trace / 3.0_kr

      DeviatoricPart3DS%YZ = M1%YZ
      DeviatoricPart3DS%XZ = M1%XZ
      DeviatoricPart3DS%XY = M1%XY
   end function DeviatoricPart3DS

   function HydrostaticPart2D(M1)
      type(Mat2D), intent(IN)                      :: M1
      type(Mat2D)                                  :: HydrostaticPart2D

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      HydrostaticPart2D%XX = M1_Trace * 0.5_kr
      HydrostaticPart2D%XY = 0.0_kr
      HydrostaticPart2D%YX = 0.0_kr
      HydrostaticPart2D%YY = M1_Trace * 0.5_kr
   end function HydrostaticPart2D

   function HydrostaticPart2DS(M1)
      type(MatS2D), intent(IN)                     :: M1
      type(MatS2D)                                 :: HydrostaticPart2DS

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      HydrostaticPart2DS%XX = M1_Trace * 0.5_kr
      HydrostaticPart2DS%YY = M1_Trace * 0.5_kr
      HydrostaticPart2DS%XY = 0.0_kr
   end function HydrostaticPart2DS

   function HydrostaticPart3D(M1)
      type(Mat3D), intent(IN)                      :: M1
      type(Mat3D)                                  :: HydrostaticPart3D

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      HydrostaticPart3D%XX = M1_Trace / 3.0_kr
      HydrostaticPart3D%XY = 0.0_kr
      HydrostaticPart3D%XZ = 0.0_kr

      HydrostaticPart3D%YX = 0.0_kr
      HydrostaticPart3D%YY = M1_Trace / 3.0_kr
      HydrostaticPart3D%YZ = 0.0_kr

      HydrostaticPart3D%ZX = 0.0_kr
      HydrostaticPart3D%ZY = 0.0_kr
      HydrostaticPart3D%ZZ = M1_Trace / 3.0_kr
   end function HydrostaticPart3D

   function HydrostaticPart3DS(M1)
      type(MatS3D), intent(IN)                     :: M1
      type(MatS3D)                                 :: HydrostaticPart3DS

      PetscReal                                    :: M1_Trace

      M1_Trace = Trace(M1)
      HydrostaticPart3DS%XX = M1_Trace / 3.0_kr
      HydrostaticPart3DS%YY = M1_Trace / 3.0_kr
      HydrostaticPart3DS%ZZ = M1_Trace / 3.0_kr

      HydrostaticPart3DS%YZ = 0.0_kr
      HydrostaticPart3DS%XZ = 0.0_kr
      HydrostaticPart3DS%XY = 0.0_kr
   end function HydrostaticPart3DS

!====================================================================
!             END OF OPERATOR OVERLOADING
!====================================================================

   function Vol_Tetra_3D(V1, V2, V3, V4)
      type(Vect3D), intent(IN)                     :: V1, V2, V3, V4
      PetscReal                                    :: Vol_Tetra_3D

      type(Vect3D)                                 :: C1, C2, C3

      C1 = V1 - V4
      C2 = V2 - V4
      C3 = V3 - V4
      Vol_Tetra_3D = abs(C1%X * (C2%Y * C3%Z - C2%Z * C3%Y) - C1%Y * (C2%X * C3%Z - C2%Z * C3%X) &
                         + C1%Z * (C2%X * C3%Y - C2%Y * C3%X)) / 6.0_kr
   end function Vol_Tetra_3D

   function Area_Tri_2D(S1, S2, S3)
      type(Vect2D), intent(IN)                     :: S1, S2, S3
      PetscReal                                    :: Area_Tri_2D

      type(Vect2D)                                 :: C1, C2

      C1 = S2 - S1
      C2 = S3 - S1
      Area_Tri_2D = abs(C1%X * C2%Y - C1%Y * C2%X) * 0.5_kr
   end function Area_Tri_2D

   function Ht_Min_Tri_2D(S1, S2, S3)
      type(Vect2D), intent(IN)                    :: S1, S2, S3
      PetscReal                                   :: Ht_Min_Tri_2D

      type(Vect2D)                                :: C1, C2, C3
      PetscReal                                   :: H1, H2, H3, AreaX2

      C1 = S2 - S3
      C2 = S3 - S1
      C3 = S1 - S2

      AreaX2 = abs(C1%X * C2%Y - C1%Y * C2%X)

      H1 = AreaX2 / sqrt((C1.DotP.C1))
      H2 = AreaX2 / sqrt((C2.DotP.C2))
      H3 = AreaX2 / sqrt((C3.DotP.C3))

      Ht_Min_Tri_2D = min(H1, H2, H3)
   end function Ht_Min_Tri_2D

   function DetMat2D(M)
      type(Mat2D), intent(IN)                     :: M
      PetscReal                                   :: DetMat2D

      DetMat2D = M%XX * M%YY - M%XY * M%YX
   end function DetMat2D

   function DetMatS2D(M)
      type(MatS2D), intent(IN)                    :: M
      PetscReal                                   :: DetMatS2D

      DetMatS2D = M%XX * M%YY - M%XY * M%XY
   end function DetMatS2D

   function DetMat3D(M)
      type(Mat3D), intent(IN)                     :: M
      PetscReal                                   :: DetMat3D

      DetMat3D = M%XX * (M%YY * M%ZZ - M%ZY * M%YZ) - M%YX * (M%XY * M%ZZ - M%ZY * M%XZ) + M%ZX * (M%XY * M%YZ - M%YY * M%XZ)
   end function DetMat3D

   function DetMatS3D(M)
      type(MatS3D), intent(IN)                    :: M
      PetscReal                                   :: DetMatS3D

      DetMatS3D = M%XX * (M%YY * M%ZZ - M%YZ * M%YZ) - M%XY * (M%XY * M%ZZ - M%YZ * M%XZ) + M%XZ * (M%XY * M%YZ - M%YY * M%XZ)
   end function DetMatS3D

   function InvertMat2D(M)
      type(Mat2D), intent(IN)                     :: M
      type(Mat2D)                                 :: InvertMat2D

      type(Mat2D)                                 :: CofMt
      PetscReal                                   :: DetM

      DetM = M%XX * M%YY - M%XY * M%YX
      CofMt%XX = M%YY
      CofMt%XY = -M%XY
      CofMt%YX = -M%YX
      CofMt%YY = M%XX

      InvertMat2D = CofMt / DetM
   end function InvertMat2D

   function InvertMatS2D(M)
      type(MatS2D), intent(IN)                    :: M
      type(MatS2D)                                :: InvertMatS2D

      type(MatS2D)                                :: CofMt
      PetscReal                                   :: DetM

      DetM = M%XX * M%YY - M%XY**2
      CofMt%XX = M%YY
      CofMt%XY = -M%XY
      CofMt%YY = M%XX

      InvertMatS2D = CofMt / DetM
   end function InvertMatS2D

   function InvertMat3D(M)
      type(Mat3D), intent(IN)                     :: M
      type(Mat3D)                                 :: InvertMat3D

      type(Mat3D)                                 :: CofMt
      PetscReal                                   :: DetM

      DetM = M%XX * (M%YY * M%ZZ - M%ZY * M%YZ) - M%YX * (M%XY * M%ZZ - M%ZY * M%XZ) + M%ZX * (M%XY * M%YZ - M%YY * M%XZ)

      CofMt%XX = M%YY * M%ZZ - M%ZY * M%YZ
      CofMt%YX = -(M%YX * M%ZZ - M%ZX * M%YZ)
      CofMt%ZX = M%YX * M%ZY - M%ZX * M%YY
      CofMt%XY = -(M%XY * M%ZZ - M%ZY * M%XZ)
      CofMt%YY = M%XX * M%ZZ - M%ZX * M%XZ
      CofMt%ZY = -(M%XX * M%ZY - M%ZX * M%XY)
      CofMt%XZ = M%XY * M%YZ - M%YY * M%XZ
      CofMt%YZ = -(M%XX * M%YZ - M%YX * M%XZ)
      CofMt%ZZ = M%XX * M%YY - M%YX * M%XY

      InvertMat3D = CofMt / DetM
   end function InvertMat3D

   function InvertMatS3D(M)
      type(MatS3D), intent(IN)                    :: M
      type(MatS3D)                                :: InvertMatS3D

      type(MatS3D)                                :: CofMt
      PetscReal                                   :: DetM

      DetM = M%XX * (M%YY * M%ZZ - M%YZ * M%YZ) - M%XY * (M%XY * M%ZZ - M%YZ * M%XZ) + M%XZ * (M%XY * M%YZ - M%YY * M%XZ)

      CofMt%XX = M%YY * M%ZZ - M%YZ * M%YZ
      CofMt%XY = -(M%XY * M%ZZ - M%XZ * M%YZ)
      CofMt%XZ = M%XY * M%YZ - M%XZ * M%YY
      CofMt%YY = M%XX * M%ZZ - M%XZ * M%XZ
      CofMt%YZ = -(M%XX * M%YZ - M%XY * M%XZ)
      CofMt%ZZ = M%XX * M%YY - M%XY * M%XY

      InvertMatS3D = CofMt / DetM
   end function InvertMatS3D

#undef __FUNCT__
#define __FUNCT__ "simplexNormal2D"
   !!! author: Blaise Bourdin (2013, bourdin@lsu.edu)
   !!!
   !!!  simplexNormal2D: Compute the normal to a simplex in 2D
   !!!
   subroutine simplexNormal2D(Coord, n, ierr)
      type(Vect2D), dimension(:), pointer              :: Coord
      type(Vect2D), intent(OUT)                        :: n
      PetscErrorCode, intent(INOUT)                    :: ierr

      n = [Coord(1)%Y - Coord(2)%Y, Coord(2)%X - Coord(1)%X]
      n = n / norm(n)
   end subroutine simplexNormal2D

#undef __FUNCT__
#define __FUNCT__ "simplexNormal3D"
   !!! author: Blaise Bourdin (2013, bourdin@lsu.edu)
   !!!
   !!!  simplexNormal3D: Compute the normal to a simplex in 3D
   !!!
   subroutine simplexNormal3D(Coord, n, ierr)
      type(Vect3D), dimension(:), pointer              :: Coord
      type(Vect3D), intent(OUT)                        :: n
      PetscErrorCode, intent(INOUT)                    :: ierr

      n = (Coord(2) - Coord(1)) .crossP. (Coord(1) - Coord(3))
      n = n / norm(n)
      ierr = 0
   end subroutine simplexNormal3D

   function InvertTens4OS2D(T)
      type(Tens4OS2D), intent(IN)                  :: T
      type(Tens4OS2D)                              :: InvertTens4OS2D

      integer                                      :: ierr

      PetscReal, dimension(3, 3)                   :: TmpArray
      PetscInt, dimension(3)                       :: ipiv
      PetscReal, dimension(3)                      :: work

      !! We convert T in a matrix using Mandel notations,invert the matrix then write back in a tensor
      TmpArray = T
      call DGETRF(3, 3, TmpArray, 3, ipiv, ierr)
      call DGETRI(3, TmpArray, 3, ipiv, work, 3, ierr)
      InvertTens4OS2D = TmpArray
   end function InvertTens4OS2D

   function InvertTens4OS3D(T)
      type(Tens4OS3D), intent(IN)                  :: T
      type(Tens4OS3D)                              :: InvertTens4OS3D

      integer                                      :: ierr

      PetscReal, dimension(6, 6)                   :: TmpArray
      PetscInt, dimension(6)                       :: ipiv
      PetscReal, dimension(6)                      :: work

      !! We convert T in a matrix using Mandel notations,invert the matrix then write back in a tensor
      TmpArray = T
      call DGETRF(6, 6, TmpArray, 6, ipiv, ierr)
      call DGETRI(6, TmpArray, 6, ipiv, work, 6, ierr)
      InvertTens4OS3D = TmpArray
   end function InvertTens4OS3D

   subroutine Tens4OS2D2Array4(A, T)
      PetscReal, dimension(2, 2, 2, 2), intent(OUT)   :: A
      type(Tens4OS2D), intent(IN)                     :: T

      A(1, 1, 1, 1) = T%XXXX
      A(1, 1, 1, 2) = T%XXXY
      A(1, 1, 2, 1) = T%XXXY
      A(1, 1, 2, 2) = T%XXYY
      A(1, 2, 1, 1) = T%XXXY
      A(1, 2, 1, 2) = T%XYXY
      A(1, 2, 2, 1) = T%XYXY
      A(1, 2, 2, 2) = T%YYXY
      A(2, 1, 1, 1) = T%XXXY
      A(2, 1, 1, 2) = T%XYXY
      A(2, 1, 2, 1) = T%XYXY
      A(2, 1, 2, 2) = T%YYXY
      A(2, 2, 1, 1) = T%XXYY
      A(2, 2, 1, 2) = T%YYXY
      A(2, 2, 2, 1) = T%YYXY
      A(2, 2, 2, 2) = T%YYYY
   end subroutine Tens4OS2D2Array4

   subroutine Tens4OS3D2Array4(A, T)
      PetscReal, dimension(3, 3, 3, 3), intent(OUT)   :: A
      type(Tens4OS3D), intent(IN)                     :: T

      A(1, 1, 1, 1) = T%XXXX
      A(1, 1, 2, 2) = T%XXYY; A(2, 2, 1, 1) = T%XXYY
      A(1, 1, 3, 3) = T%XXZZ; A(3, 3, 1, 1) = T%XXZZ
      A(1, 1, 2, 3) = T%XXYZ; A(1, 1, 3, 2) = T%XXYZ; A(2, 3, 1, 1) = T%XXYZ; A(3, 2, 1, 1) = T%XXYZ
      A(1, 1, 1, 3) = T%XXXZ; A(1, 1, 3, 1) = T%XXXZ; A(1, 3, 1, 1) = T%XXXZ; A(3, 1, 1, 1) = T%XXXZ
      A(1, 1, 1, 2) = T%XXXY; A(1, 1, 2, 1) = T%XXXY; A(1, 2, 1, 1) = T%XXXY; A(2, 1, 1, 1) = T%XXXY

      A(2, 2, 2, 2) = T%YYYY
      A(2, 2, 3, 3) = T%YYZZ; A(3, 3, 2, 2) = T%YYZZ
      A(2, 2, 2, 3) = T%YYYZ; A(2, 2, 3, 2) = T%YYYZ; A(2, 3, 2, 2) = T%YYYZ; A(3, 2, 2, 2) = T%YYYZ
      A(2, 2, 1, 3) = T%YYXZ; A(2, 2, 3, 1) = T%YYXZ; A(1, 3, 2, 2) = T%YYXZ; A(3, 1, 2, 2) = T%YYXZ
      A(2, 2, 1, 2) = T%YYXY; A(2, 2, 2, 1) = T%YYXY; A(1, 2, 2, 2) = T%YYXY; A(2, 1, 2, 2) = T%YYXY

      A(3, 3, 3, 3) = T%ZZZZ
      A(3, 3, 2, 3) = T%ZZYZ; A(3, 3, 3, 2) = T%ZZYZ; A(2, 3, 3, 3) = T%ZZYZ; A(3, 2, 3, 3) = T%ZZYZ
      A(3, 3, 1, 3) = T%ZZXZ; A(3, 3, 3, 1) = T%ZZXZ; A(1, 3, 3, 3) = T%ZZXZ; A(3, 1, 3, 3) = T%ZZXZ
      A(3, 3, 1, 2) = T%ZZXY; A(3, 3, 2, 1) = T%ZZXY; A(1, 2, 3, 3) = T%ZZXY; A(2, 1, 3, 3) = T%ZZXY

      A(2, 3, 2, 3) = T%YZYZ; A(2, 3, 3, 2) = T%YZYZ; A(3, 2, 2, 3) = T%YZYZ; A(3, 2, 3, 2) = T%YZYZ
      A(2, 3, 1, 3) = T%YZXZ; A(2, 3, 3, 1) = T%YZXZ; A(3, 2, 1, 3) = T%YZXZ; A(3, 2, 3, 1) = T%YZXZ
      A(1, 3, 2, 3) = T%YZXZ; A(1, 3, 3, 2) = T%YZXZ; A(3, 1, 2, 3) = T%YZXZ; A(3, 1, 3, 2) = T%YZXZ
      A(2, 3, 1, 2) = T%YZXY; A(2, 3, 2, 1) = T%YZXY; A(3, 2, 1, 2) = T%YZXY; A(3, 2, 2, 1) = T%YZXY
      A(1, 2, 2, 3) = T%YZXY; A(1, 2, 3, 2) = T%YZXY; A(2, 1, 2, 3) = T%YZXY; A(2, 1, 3, 2) = T%YZXY

      A(1, 2, 1, 2) = T%XYXY; A(1, 2, 2, 1) = T%XYXY; A(2, 1, 1, 2) = T%XYXY; A(2, 1, 2, 1) = T%XYXY
      A(1, 3, 1, 2) = T%XZXY; A(1, 3, 2, 1) = T%XZXY; A(3, 1, 1, 2) = T%XZXY; A(3, 1, 2, 1) = T%XZXY
      A(1, 2, 1, 3) = T%XZXY; A(1, 2, 3, 1) = T%XZXY; A(2, 1, 1, 3) = T%XZXY; A(2, 1, 3, 1) = T%XZXY
      A(1, 3, 1, 3) = T%XZXZ; A(1, 3, 3, 1) = T%XZXZ; A(3, 1, 1, 3) = T%XZXZ; A(3, 1, 3, 1) = T%XZXZ
   end subroutine Tens4OS3D2Array4

   subroutine Array42Tens4OS2D(T, A)
      type(Tens4OS2D), intent(OUT)                    :: T
      PetscReal, dimension(2, 2, 2, 2), intent(IN)    :: A

      T%XXXX = A(1, 1, 1, 1)
      T%XXXY = A(1, 1, 1, 2)
      T%XXYY = A(1, 1, 2, 2)
      T%XYXY = A(1, 2, 1, 2)
      T%YYXY = A(2, 2, 1, 2)
      T%YYYY = A(2, 2, 2, 2)
   end subroutine Array42Tens4OS2D

   subroutine Array42Tens4OS3D(T, A)
      type(Tens4OS3D), intent(OUT)                    :: T
      PetscReal, dimension(3, 3, 3, 3), intent(IN)    :: A

      T%XXXX = A(1, 1, 1, 1); T%XXYY = A(1, 1, 2, 2); T%XXZZ = A(1, 1, 3, 3); T%XXYZ = A(1, 1, 2, 3); T%XXXZ = A(1, 1, 1, 3); T%XXXY = A(1, 1, 1, 2)
      T%YYYY = A(2, 2, 2, 2); T%YYZZ = A(2, 2, 3, 3); T%YYYZ = A(2, 2, 2, 3); T%YYXZ = A(2, 2, 1, 3); T%YYXY = A(2, 2, 1, 2)
      T%ZZZZ = A(3, 3, 3, 3); T%ZZYZ = A(3, 3, 2, 3); T%ZZXZ = A(3, 3, 1, 3); T%ZZXY = A(3, 3, 1, 2)
      T%YZYZ = A(2, 3, 2, 3); T%YZXZ = A(2, 3, 1, 3); T%YZXY = A(2, 3, 1, 2)
      T%XZXZ = A(1, 3, 1, 3); T%XZXY = A(1, 3, 1, 2)
      T%XYXY = A(1, 2, 1, 2)
   end subroutine Array42Tens4OS3D

   function Tens4OS2DTransform(T, M)
      !!! Apply the transformation given by the matrix R to a 4th order tensor
      !!! i.e. C_{ijkl} = R_{ip}.R_{jq}.R_{kr}.R{ls} A_{pqrs}
      type(Tens4OS2D), intent(IN)                  :: T
      type(Mat2D), intent(IN)                      :: M

      type(Tens4OS2D)                              :: Tens4OS2DTransform

      PetscReal, dimension(2, 2, 2, 2)             :: TT, C
      PetscReal, dimension(2, 2)                   :: MM
      integer                                      :: i, j, k, l
      integer                                      :: p, q, r, s

      TT = T
      MM = M
      C = 0.0_kr
      do i = 1, 2
         do j = 1, 2
            do k = 1, 2
               do l = 1, 2
                  do p = 1, 2
                     do q = 1, 2
                        do r = 1, 2
                           do s = 1, 2
                              C(i, j, k, l) = C(i, j, k, l) + MM(i, p) * MM(j, q) * MM(k, r) * MM(l, s) * TT(p, q, r, s)
                           end do
                        end do
                     end do
                  end do
               end do
            end do
         end do
      end do

      Tens4OS2DTransform = C

   end function Tens4OS2DTransform

   function Tens4OS3DTransform(T, M)
      !!! Apply the transformation given by the matrix R to a 4th order tensor
      !!! i.e. C_{ijkl} = R_{ip}.R_{jq}.R_{kr}.R{ls} A_{pqrs}
      type(Tens4OS3D), intent(IN)                 :: T
      type(Mat3D), intent(IN)                     :: M

      type(Tens4OS3D)                             :: Tens4OS3DTransform

      PetscReal, dimension(3, 3, 3, 3)            :: TT, C
      PetscReal, dimension(3, 3)                  :: MM
      integer                                     :: i, j, k, l
      integer                                     :: p, q, r, s

      TT = T
      MM = M
      C = 0.0_kr
      do i = 1, 3
         do j = 1, 3
            do k = 1, 3
               do l = 1, 3
                  do p = 1, 3
                     do q = 1, 3
                        do r = 1, 3
                           do s = 1, 3
                              C(i, j, k, l) = C(i, j, k, l) + MM(i, p) * MM(j, q) * MM(k, r) * MM(l, s) * TT(p, q, r, s)
                           end do
                        end do
                     end do
                  end do
               end do
            end do
         end do
      end do

      Tens4OS3DTransform = C
   end function Tens4OS3DTransform

   function Tens4OS2DSquareRoot(T)
      type(Tens4OS2D), intent(IN)                  :: T
      type(Tens4OS2D)                              :: Tens4OS2DSquareRoot

      integer, parameter                           :: n = 3
      integer                                      :: i, j
      PetscReal, dimension(n, n)                   :: A, Pt
      PetscReal, dimension(n)                      :: lmbda
      PetscReal                                    :: d
      PetscInt                                     :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)       :: work
      PetscInt                                     :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                 :: iwork
      PetscInt                                     :: info

      A = T
      call DSYEVD('V', 'L', n, A, n, lmbda, work, lwork, iwork, liwork, info)
      Pt = transpose(A)
      do i = 1, n

         if (lmbda(i) < 0.0_kr) then
            write (*, *) 'ERROR in Tens4OSDSquareRoot, negative eigenvalue ', lmbda(i)
         else
            d = sqrt(lmbda(i))
            do j = 1, n
               Pt(i, j) = A(j, i) * d
            end do
         end if
      end do
      Tens4OS2DSquareRoot = matmul(A, Pt)
   end function Tens4OS2DSquareRoot

   function Tens4OS3DSquareRoot(T)
      type(Tens4OS3D), intent(IN)                  :: T
      type(Tens4OS3D)                              :: Tens4OS3DSquareRoot

      integer, parameter                           :: n = 6
      integer                                      :: i, j
      PetscReal, dimension(n, n)                   :: A, Pt
      PetscReal, dimension(n)                      :: lmbda
      PetscReal                                    :: d
      PetscInt                                     :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)       :: work
      PetscInt                                     :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                 :: iwork
      PetscInt                                     :: info

      A = T
      call DSYEVD('V', 'L', n, A, n, lmbda, work, lwork, iwork, liwork, info)
      Pt = transpose(A)
      do i = 1, n
         if (lmbda(i) < 0.0_kr) then
            write (*, *) 'ERROR in Tens4OSDSquareRoot, negative eigenvalue ', lmbda(i)
         else
            d = sqrt(lmbda(i))
            do j = 1, n
               Pt(i, j) = A(j, i) * d
            end do
         end if
      end do
      Tens4OS3DSquareRoot = matmul(A, Pt)
   end function Tens4OS3DSquareRoot

   subroutine MatS2DSpectralDecomposition(M, ppleValues, ppleDirections)
      type(MatS2D), intent(IN)                      :: M
      PetscReal, dimension(2), intent(OUT)          :: ppleValues
      type(MatS2D), dimension(2), intent(OUT)       :: ppleDirections

      integer, parameter                            :: n = 2
      PetscReal, dimension(n, n)                    :: A
      PetscInt                                      :: i
      PetscInt                                      :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)        :: work
      PetscInt                                      :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                  :: iwork
      PetscInt                                      :: info

      A = M
      call DSYEVD('V', 'L', n, A, n, ppleValues, work, lwork, iwork, liwork, info)
      do i = 1, n
         ppleDirections(i)%XX = A(1, i)**2
         ppleDirections(i)%YY = A(2, i)**2
         ppleDirections(i)%XY = A(1, i) * A(2, i)
      end do
   end subroutine MatS2DSpectralDecomposition

   subroutine MatS3DSpectralDecomposition(M, ppleValues, ppleDirections)
      type(MatS3D), intent(IN)                      :: M
      PetscReal, dimension(3), intent(OUT)          :: ppleValues
      type(MatS3D), dimension(3), intent(OUT)       :: ppleDirections

      integer, parameter                            :: n = 3
      PetscReal, dimension(n, n)                    :: A
      PetscInt                                      :: i
      PetscInt                                      :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)        :: work
      PetscInt                                      :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                  :: iwork
      PetscInt                                      :: info

      A = M
      call DSYEVD('V', 'L', n, A, n, ppleValues, work, lwork, iwork, liwork, info)
      do i = 1, n
         ppleDirections(i)%XX = A(1, i)**2
         ppleDirections(i)%YY = A(2, i)**2
         ppleDirections(i)%ZZ = A(3, i)**2
         ppleDirections(i)%YZ = A(2, i) * A(3, i)
         ppleDirections(i)%XZ = A(1, i) * A(3, i)
         ppleDirections(i)%XY = A(1, i) * A(2, i)
      end do
   end subroutine MatS3DSpectralDecomposition

   subroutine MatS3DEigenVectorValues(M, MatProj, MatDiag)
      type(MatS3D), intent(IN)                     :: M
      PetscReal, dimension(3)                      :: ppleValues
      type(Mat3D), intent(OUT)                     :: MatProj
      type(MatS3D), intent(OUT)                    :: MatDiag

      integer, parameter                           :: n = 3
      PetscReal, dimension(n, n)                   :: A
      PetscInt                                     :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)       :: work
      PetscInt                                     :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                 :: iwork
      PetscInt                                     :: info

      A = M
      call DSYEVD('V', 'L', n, A, n, ppleValues, work, lwork, iwork, liwork, info)

      MatDiag = 0.0_kr
      MatDiag%XX = ppleValues(1)
      MatDiag%YY = ppleValues(2)
      MatDiag%ZZ = ppleValues(3)

      MatProj%XX = A(1, 1)
      MatProj%YX = A(2, 1)
      MatProj%ZX = A(3, 1)
      MatProj%XY = A(1, 2)
      MatProj%YY = A(2, 2)
      MatProj%ZY = A(3, 2)
      MatProj%XZ = A(1, 3)
      MatProj%YZ = A(2, 3)
      MatProj%ZZ = A(3, 3)
   end subroutine MatS3DEigenVectorValues

   subroutine MatS2DEigenVectorValues(M, MatProj, MatDiag)
      type(MatS2D), intent(IN)                     :: M
      PetscReal, dimension(2)                      :: ppleValues
      type(Mat2D), intent(OUT)                     :: MatProj
      type(MatS2D), intent(OUT)                    :: MatDiag

      integer, parameter                           :: n = 2
      PetscReal, dimension(n, n)                   :: A
      PetscInt                                     :: lwork = 2 * n**2 + 6 * n + 1
      PetscReal, dimension(2*n**2 + 6*n + 1)       :: work
      PetscInt                                     :: liwork = 5 * n + 3
      PetscInt, dimension(5*n + 3)                 :: iwork
      PetscInt                                     :: info

      A = M
      call DSYEVD('V', 'L', n, A, n, ppleValues, work, lwork, iwork, liwork, info)
      if (info /= 0) then
         write (*, *) 'DSYEVD failed with info=', info
         write (*, *) 'A: ', A
         write (*, *) 'ppleValues: ', ppleValues
         SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_LIB, "DSYEVD failed: "//__FUNCT__)
      end if

      MatDiag = 0.0_kr
      MatDiag%XX = ppleValues(1)
      MatDiag%YY = ppleValues(2)

      MatProj%XX = A(1, 1)
      MatProj%XY = A(1, 2)
      MatProj%YX = A(2, 1)
      MatProj%YY = A(2, 2)
   end subroutine MatS2DEigenVectorValues

#undef __FUNCT__
#define __FUNCT__ "Mat2DMoment"
!!! author: Blaise Bourdin (2019, bourdin@lsu.edu)
!!!
!!!  Mat2DMoment: k-th moment of a 2x2 matrix, i.e tr(A^k)
!!!

   function Mat2DMoment(k, A)
      integer, intent(IN)                          :: k
      type(Mat2D), intent(IN)                      :: A
      PetscReal                                    :: Mat2DMoment

      integer                                      :: i
      type(Mat2D)                                  :: Ak

      select case (k)
      case (1)
         Mat2DMoment = trace(A)
      case (2)
         Mat2DMoment = (A%XX**2 + 2 * A%XY * A%YX + A%YY**2) / 2.0_kr
      case default
         Ak = A
         do i = 1, k - 1
            Ak = Ak * A
         end do
         Mat2DMoment = trace(Ak) / k
      end select
   end function Mat2DMoment

#undef __FUNCT__
#define __FUNCT__ "MatS2DMoment"
!!! author: Blaise Bourdin (2019, bourdin@lsu.edu)
!!!
!!!  MatS2DMoment: k-th moment of a 2x2 symmetric matrix, i.e tr(A^k)
!!!

   function MatS2DMoment(k, A)
      integer, intent(IN)                         :: k
      type(MatS2D), intent(IN)                    :: A
      PetscReal                                   :: MatS2DMoment

      integer                                     :: i
      type(MatS2D)                                :: Ak

      select case (k)
      case (1)
         MatS2DMoment = trace(A)
      case (2)
         MatS2DMoment = (A%XX**2 + 2 * A%XY**2 + A%YY**2) / 2.0_kr
      case default
         Ak = A
         do i = 1, k - 1
            Ak = Ak * A
         end do
         MatS2DMoment = trace(Ak) / k
      end select
   end function MatS2DMoment

#undef __FUNCT__
#define __FUNCT__ "Mat3DMoment"
!!! author: Blaise Bourdin (2019, bourdin@lsu.edu)
!!!
!!!  Mat3DMoment: k-th moment of a 3x3 matrix, i.e tr(A^k)
!!!

   function Mat3DMoment(k, A)
      integer, intent(IN)                         :: k
      type(Mat3D), intent(IN)                     :: A
      PetscReal                                   :: Mat3DMoment

      integer                                     :: i
      type(Mat3D)                                 :: Ak

      select case (k)
      case (1)
         Mat3DMoment = trace(A)
      case (2)
         Mat3DMoment = (A%XX**2 + A%YY**2 + A%ZZ**2 + 2 * A%XY * A%YX + 2 * A%XZ * A%ZX + 2 * A%YZ * A%ZY) / 2.0_kr
      case default
         Ak = A
         do i = 1, k - 1
            Ak = Ak * A
         end do
         Mat3DMoment = trace(Ak) / k
      end select
   end function Mat3DMoment

#undef __FUNCT__
#define __FUNCT__ "MatS3DMoment"
!!! author: Blaise Bourdin (2019, bourdin@lsu.edu)
!!!
!!!  MatS3DMoment: k-th moment of a 3x3 symmetric matrix, i.e tr(A^k)
!!!

   function MatS3DMoment(k, A)
      integer, intent(IN)                         :: k
      type(MatS3D), intent(IN)                    :: A
      PetscReal                                   :: MatS3DMoment

      integer                                     :: i
      type(MatS3D)                                :: Ak

      select case (k)
      case (1)
         MatS3DMoment = trace(A)
      case (2)
         MatS3DMoment = (A%XX**2 + A%YY**2 + A%ZZ**2 + 2 * A%YZ**2 + 2 * A%XZ**2 + 2 * A%XY**2) / 2.0_kr
      case default
         Ak = A
         do i = 1, k - 1
            Ak = Ak * A
         end do
         MatS3DMoment = trace(Ak) / k
      end select
   end function MatS3DMoment
end module m_MEF90_LinAlg
