1!/*****************************************************************************/ 2! * 3! * Elmer/Ice, a glaciological add-on to Elmer 4! * http://elmerice.elmerfem.org 5! * 6! * 7! * This program is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU General Public License 9! * as published by the Free Software Foundation; either version 2 10! * of the License, or (at your option) any later version. 11! * 12! * This program is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! * GNU General Public License for more details. 16! * 17! * You should have received a copy of the GNU General Public License 18! * along with this program (in file fem/GPL-2); if not, write to the 19! * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 20! * Boston, MA 02110-1301, USA. 21! * 22! *****************************************************************************/ 23! ****************************************************************************** 24! * 25! * Authors: Fabien / OG 26! * Email: 27! * Web: http://elmerice.elmerfem.org 28! * 29! * Original Date: 13/10/05 from version 1.5 30! * 31! *****************************************************************************/ 32!------------------------------------------------------------------------------ 33 RECURSIVE SUBROUTINE ComputeEigenValues( Model,Solver,dt,TransientSimulation ) 34!------------------------------------------------------------------------------ 35 36 USE DefUtils 37 38 IMPLICIT NONE 39 40!------------------------------------------------------------------------------ 41!****************************************************************************** 42! 43! ARGUMENTS: 44! 45! TYPE(Model_t) :: Model, 46! INPUT: All model information (mesh,materials,BCs,etc...) 47! 48! TYPE(Solver_t) :: Solver 49! INPUT: Linear equation solver options 50! 51! REAL(KIND=dp) :: dt, 52! INPUT: Timestep size for time dependent simulations (NOTE: Not used 53! currently) 54! 55!****************************************************************************** 56 57 TYPE(Model_t) :: Model 58 TYPE(Solver_t), TARGET :: Solver 59 60 LOGICAL :: TransientSimulation 61 REAL(KIND=dp) :: dt 62!------------------------------------------------------------------------------ 63! Local variables 64!------------------------------------------------------------------------------ 65 66 67 68 INTEGER :: dim, StressDOFs 69 TYPE(Variable_t), POINTER :: StressSol,EigenSol,EigenV1,EigenV2,EigenV3 70 REAL(KIND=dp), POINTER :: Stress(:),Eigen(:) 71 INTEGER, POINTER :: StressPerm(:),EigenPerm(:) 72 73 REAL(KIND=dp),dimension(3,3) :: localM,EigenVec 74 REAL(KIND=dp),dimension(3) :: EigValues, EI, ki 75 REAL(KIND=dp) :: WORK(24),Dumy(1) 76 Real(KIND=dp) :: a 77 INTEGER :: i, j, t, ordre(3),infor 78 LOGICAL :: GotIt,UnFoundFatal=.TRUE. 79 CHARACTER(LEN=MAX_NAME_LEN) :: TensorVarName, EigenVarName 80 81 82!------------------------------------------------------------------------------ 83! Get variables needed for solution 84!------------------------------------------------------------------------------ 85 dim = CoordinateSystemDimension() 86 87 TensorVarName = GetString( Solver % Values, 'Tensor Variable Name', GotIt ) 88 IF (.NOT.Gotit) TensorVarName = 'Stress' 89 StressSol => VariableGet( Solver % Mesh % Variables, TensorVarName ) 90 StressPerm => StressSol % Perm 91 StressDOFs = StressSol % DOFs 92 Stress => StressSol % Values 93 94 ! Eigen Values (dimension 3) 95 EigenVarName = GetString( Solver % Values, 'EigenValue Variable Name', GotIt ) 96 IF (.NOT.Gotit) EigenVarName = 'EigenStress' 97 EigenSol => VariableGet( Solver % Mesh % Variables, EigenVarName,UnFoundFatal=UnFoundFatal) 98 EigenPerm => EigenSol % Perm 99 Eigen => EigenSol % Values 100 101 ! Eigen Vector (dimension 3) 102 EigenV1 => VariableGet( Solver % Mesh % Variables, 'EigenVector1' ) 103 EigenV2 => VariableGet( Solver % Mesh % Variables, 'EigenVector2' ) 104 EigenV3 => VariableGet( Solver % Mesh % Variables, 'EigenVector3' ) 105 106 107 Do t=1,Solver % Mesh % NumberOfNodes 108 109 ! the Stress components [Sxx, Syy, Szz, Sxy, Syz, Szx] 110 localM=0.0_dp 111 localM(1,1)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 1) 112 localM(2,2)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 2) 113 localM(3,3)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 3) 114 localM(1,2)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 4) 115 localM(2,1)=localM(1,2) 116 if (dim.eq.3) then 117 localM(2,3)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 5) 118 localM(1,3)=Stress( StressDOFs * (StressPerm(t) - 1 ) + 6) 119 localM(3,2)=localM(2,3) 120 localM(3,1)=localM(1,3) 121 end if 122 123! Compute EigenValues using lapack DGEEV subroutine 124 CALL DGEEV('N','V',3,localM,3,EigValues,EI,Dumy,1,EigenVec,3,Work,24,infor ) 125 IF (infor.ne.0) & 126 CALL FATAL('Compute EigenValues', 'Failed to compute EigenValues') 127 localM(1:3,1:3)=EigenVec(1:3,1:3) 128! Ordered value Ev1 < Ev2 < Ev3 129 130 Do i=1,3 131 ki(i)=EigValues(i) 132 ordre(i)=i 133 End Do 134 Do j=2,3 135 a=ki(j) 136 Do i=j-1,1,-1 137 If (ki(i).LE.a) Goto 20 138 ki(i+1)=ki(i) 139 ordre(i+1)=ordre(i) 140 End Do 141 20 Continue 142 ki(i+1)=a 143 ordre(i+1)=j 144 End Do 145 146 Eigen( 3 * ( EigenPerm(t) - 1) + 1) = EigValues(ordre(1)) 147 Eigen( 3 * ( EigenPerm(t) - 1) + 2) = EigValues(ordre(2)) 148 Eigen( 3 * ( EigenPerm(t) - 1) + 3) = EigValues(ordre(3)) 149 150 If (associated(EigenV1)) then 151 EigenV1 % Values( 3 * (EigenV1 % Perm(t) - 1) + 1 ) = localM(ordre(1),1) 152 EigenV1 % Values( 3 * (EigenV1 % Perm(t) - 1) + 2 ) = localM(ordre(1),2) 153 EigenV1 % Values( 3 * (EigenV1 % Perm(t) - 1) + 3 ) = localM(ordre(1),3) 154 endif 155 If (associated(EigenV2)) then 156 EigenV2 % Values( 3 * (EigenV2 % Perm(t) - 1) + 1 ) = localM(ordre(2),1) 157 EigenV2 % Values( 3 * (EigenV2 % Perm(t) - 1) + 2 ) = localM(ordre(2),2) 158 EigenV2 % Values( 3 * (EigenV2 % Perm(t) - 1) + 3 ) = localM(ordre(2),3) 159 endif 160 If (associated(EigenV3)) then 161 EigenV3 % Values( 3 * (EigenV3 % Perm(t) - 1) + 1 ) = localM(ordre(3),1) 162 EigenV3 % Values( 3 * (EigenV3 % Perm(t) - 1) + 2 ) = localM(ordre(3),2) 163 EigenV3 % Values( 3 * (EigenV3 % Perm(t) - 1) + 3 ) = localM(ordre(3),3) 164 Endif 165 166 End do 167 168!------------------------------------------------------------------------------ 169 END SUBROUTINE ComputeEigenValues 170!------------------------------------------------------------------------------ 171 172