femwind/fortran/module_hexa.f90

   1 module module_hexa   ! testing only
   2 use module_utils
   3
   4 contains
   5
   6 subroutine hexa(A,X,u0,Kloc,Floc,Jg,iflag)
   7 ! purpose: create local stiffness matrix etc for hexa 3d element
   8 ! in:
   9 !   A   coefficient matrix size 3x3, symmetric positive definite
  10 !   X   nodes coordinates size 3x8, one each column is one node
  11 !   u0   column vector of input wind at the center of the element
  12 !   iflags  iflags = 1 compute Kloc, iflag = 2 compute Floc, iflag = 3 compute Jg
  13 ! out:
  14 !   Kloc   local stiffness matrix
  15 !   Floc   local divergence load vector
  16 !   Jg     gradient at center of function with values V is V'*Jg
  17
  18 implicit none
  19
  20 !*** arguments
  21
  22 real, intent(in):: A(3,3), X(3,8), u0(3)    ! fortran is not case sensitive
  23 integer, intent(in)::iflag
  24 real, intent(out):: Kloc(8,8), Floc(8), Jg(8,3)
  25 !*** local variables
  26 !real, parameter :: g = 0.5773502691896257
  27 real,dimension(9,3),save :: ib = reshape((/-1,-1,-1,-1,1,1,1,1,0,-1,-1,1,1,-1,-1,1,1,0,-1,1,-1,1,-1,1,-1,1,0/),(/9,3/))
  28 real,dimension(9,3),save :: s = reshape((/-0.5773502691896257,-0.5773502691896257,-0.5773502691896257,-0.5773502691896257&
  29 ,0.5773502691896257,0.5773502691896257,0.5773502691896257,0.5773502691896257,0.0,-0.5773502691896257,-0.5773502691896257&
  30 ,0.5773502691896257,0.5773502691896257,-0.5773502691896257,-0.5773502691896257,0.5773502691896257,0.5773502691896257&
  31 ,0.0,-0.5773502691896257,0.5773502691896257,-0.5773502691896257,0.5773502691896257,-0.5773502691896257,0.5773502691896257&
  32 ,-0.5773502691896257,0.5773502691896257,0.0/),(/9,3/))
  33 real :: gradf(8,3), Jg_tmp(8,3)
  34 real :: Jx(3,3), Jx_inv(3,3)
  35 real :: Jg_tran(3,8), A_tmp(3,8)
  36 real :: K_at_s(8,8)
  37 real :: tmp_mat(8)
  38 real :: u0_tmp(3)
  39 integer :: i,j,k,m
  40 real :: detJx = 0
  41 real :: tmp = 0
  42 real :: vol = 0
  43 gradf = 0.
  44 Jg_tmp = 0.
  45 Jx = 0.
  46 Jx_inv = 0.
  47 Jg_tran = 0.
  48 A_tmp = 0.
  49 K_at_s = 0.
  50 tmp_mat = 0.
  51 u0_tmp = 0.
  52
  53 Kloc = 0.
  54 Floc = 0.
  55 Jg = 0.
  56 !*** executable
  57
  58 !Calculate Jg loop
  59 if (iflag .eq.  3) then
  60
  61     do i=1,9
  62
  63     !Calculate gradf
  64     !first column of gradf
  65     do j = 1,8
  66     gradf(j,1) = (ib(j,1)*(1+ib(j,2)*s(i,2))*(1+ib(j,3)*s(i,3)))/8
  67     end do
  68     !second coumn of gradf
  69     do j = 1,8
  70     gradf(j,2) = ((1+ib(j,1)*s(i,1))*ib(j,2)*(1+ib(j,3)*s(i,3)))/8
  71     end do
  72     !third column of gradf
  73     do j = 1,8
  74     gradf(j,3) = ((1+ib(j,1)*s(i,1))*(1+ib(j,2)*s(i,2))*ib(j,3))/8
  75     end do
  76
  77     Jx = matmul(X,gradf)
  78
  79     !!!Compute Jx_inv!!!
  80     detJx = (Jx(1,1)*Jx(2,2)*Jx(3,3) - Jx(1,1)*Jx(2,3)*Jx(3,2)-&
  81              Jx(1,2)*Jx(2,1)*Jx(3,3) + Jx(1,2)*Jx(2,3)*Jx(3,1)+&
  82              Jx(1,3)*Jx(2,1)*Jx(3,2) - Jx(1,3)*Jx(2,2)*Jx(3,1))
  83
  84     if(abs(detJx).lt.tiny(detJx))then
  85         print *,'detJx=',detJx
  86         call crash('The Jacobian is (numerically) singular')
  87     endif
  88
  89     Jx_inv(1,1) = +(1/detJx) * (Jx(2,2)*Jx(3,3) - Jx(2,3)*Jx(3,2))
  90     Jx_inv(2,1) = -(1/detJx) * (Jx(2,1)*Jx(3,3) - Jx(2,3)*Jx(3,1))
  91     Jx_inv(3,1) = +(1/detJx) * (Jx(2,1)*Jx(3,2) - Jx(2,2)*Jx(3,1))
  92     Jx_inv(1,2) = -(1/detJx) * (Jx(1,2)*Jx(3,3) - Jx(1,3)*Jx(3,2))
  93     Jx_inv(2,2) = +(1/detJx) * (Jx(1,1)*Jx(3,3) - Jx(1,3)*Jx(3,1))
  94     Jx_inv(3,2) = -(1/detJx) * (Jx(1,1)*Jx(3,2) - Jx(1,2)*Jx(3,1))
  95     Jx_inv(1,3) = +(1/detJx) * (Jx(1,2)*Jx(2,3) - Jx(1,3)*Jx(2,2))
  96     Jx_inv(2,3) = -(1/detJx) * (Jx(1,1)*Jx(2,3) - Jx(1,3)*Jx(2,1))
  97     Jx_inv(3,3) = +(1/detJx) * (Jx(1,1)*Jx(2,2) - Jx(1,2)*Jx(2,1))
  98
  99     Jg = matmul(gradf,Jx_inv)
 100
 101     end do
 102 endif !end for computing Jg
 103
 104 !Calculate Kloc loop
 105 if (iflag .eq.  1) then
 106
 107 do i=1,9
 108
 109 !Calculate gradf
 110 !first column of gradf
 111 do j = 1,8
 112 gradf(j,1) = (ib(j,1)*(1+ib(j,2)*s(i,2))*(1+ib(j,3)*s(i,3)))/8
 113 end do
 114 !second coumn of gradf
 115 do j = 1,8
 116 gradf(j,2) = ((1+ib(j,1)*s(i,1))*ib(j,2)*(1+ib(j,3)*s(i,3)))/8
 117 end do
 118 !third column of gradf
 119 do j = 1,8
 120 gradf(j,3) = ((1+ib(j,1)*s(i,1))*(1+ib(j,2)*s(i,2))*ib(j,3))/8
 121 end do
 122
 123 Jx = matmul(X,gradf)
 124
 125 !!!Compute Jx_inv!!!
 126 detJx = (Jx(1,1)*Jx(2,2)*Jx(3,3) - Jx(1,1)*Jx(2,3)*Jx(3,2)-&
 127          Jx(1,2)*Jx(2,1)*Jx(3,3) + Jx(1,2)*Jx(2,3)*Jx(3,1)+&
 128          Jx(1,3)*Jx(2,1)*Jx(3,2) - Jx(1,3)*Jx(2,2)*Jx(3,1))
 129
 130     if(abs(detJx).lt.tiny(detJx))then
 131         print *,'detJx=',detJx
 132         call crash('The Jacobian is (numerically) singular')
 133     endif
 134
 135 Jx_inv(1,1) = +(1/detJx) * (Jx(2,2)*Jx(3,3) - Jx(2,3)*Jx(3,2))
 136 Jx_inv(2,1) = -(1/detJx) * (Jx(2,1)*Jx(3,3) - Jx(2,3)*Jx(3,1))
 137 Jx_inv(3,1) = +(1/detJx) * (Jx(2,1)*Jx(3,2) - Jx(2,2)*Jx(3,1))
 138 Jx_inv(1,2) = -(1/detJx) * (Jx(1,2)*Jx(3,3) - Jx(1,3)*Jx(3,2))
 139 Jx_inv(2,2) = +(1/detJx) * (Jx(1,1)*Jx(3,3) - Jx(1,3)*Jx(3,1))
 140 Jx_inv(3,2) = -(1/detJx) * (Jx(1,1)*Jx(3,2) - Jx(1,2)*Jx(3,1))
 141 Jx_inv(1,3) = +(1/detJx) * (Jx(1,2)*Jx(2,3) - Jx(1,3)*Jx(2,2))
 142 Jx_inv(2,3) = -(1/detJx) * (Jx(1,1)*Jx(2,3) - Jx(1,3)*Jx(2,1))
 143 Jx_inv(3,3) = +(1/detJx) * (Jx(1,1)*Jx(2,2) - Jx(1,2)*Jx(2,1))
 144 Jg = matmul(gradf,Jx_inv)
 145
 146 !check to calc Kloc
 147 if (i < 9) then
 148 do j = 1,8
 149 do k = 1,3
 150 Jg_tran(k,j) = Jg(j,k)
 151 end do
 152 enddo
 153
 154 A_tmp = matmul(A, Jg_tran)
 155 K_at_s = matmul(Jg,A_tmp)
 156 Kloc = Kloc+(K_at_s*abs(detJx))
 157 end if !end for computing Kloc
 158 end do !end of outter most do loop
 159 end if !end for computing Kloc
 160
 161
 162 !Calculate Floc loop
 163 if (iflag .eq. 2) then
 164
 165 do i=1,9
 166
 167 !Calculate gradf
 168 !first column of gradf
 169 do j = 1,8
 170 gradf(j,1) = (ib(j,1)*(1+ib(j,2)*s(i,2))*(1+ib(j,3)*s(i,3)))/8
 171 end do
 172 !second coumn of gradf
 173 do j = 1,8
 174 gradf(j,2) = ((1+ib(j,1)*s(i,1))*ib(j,2)*(1+ib(j,3)*s(i,3)))/8
 175 end do
 176 !third column of gradf
 177 do j = 1,8
 178 gradf(j,3) = ((1+ib(j,1)*s(i,1))*(1+ib(j,2)*s(i,2))*ib(j,3))/8
 179 end do
 180
 181 Jx = matmul(X,gradf)
 182
 183 !!!Compute Jx_inv!!!
 184 detJx = (Jx(1,1)*Jx(2,2)*Jx(3,3) - Jx(1,1)*Jx(2,3)*Jx(3,2)-&
 185          Jx(1,2)*Jx(2,1)*Jx(3,3) + Jx(1,2)*Jx(2,3)*Jx(3,1)+&
 186          Jx(1,3)*Jx(2,1)*Jx(3,2) - Jx(1,3)*Jx(2,2)*Jx(3,1))
 187
 188     if(abs(detJx).lt.tiny(detJx))then
 189         print *,'detJx=',detJx
 190         call crash('The Jacobian is (numerically) singular')
 191     endif
 192
 193 Jx_inv(1,1) = +(1/detJx) * (Jx(2,2)*Jx(3,3) - Jx(2,3)*Jx(3,2))
 194 Jx_inv(2,1) = -(1/detJx) * (Jx(2,1)*Jx(3,3) - Jx(2,3)*Jx(3,1))
 195 Jx_inv(3,1) = +(1/detJx) * (Jx(2,1)*Jx(3,2) - Jx(2,2)*Jx(3,1))
 196 Jx_inv(1,2) = -(1/detJx) * (Jx(1,2)*Jx(3,3) - Jx(1,3)*Jx(3,2))
 197 Jx_inv(2,2) = +(1/detJx) * (Jx(1,1)*Jx(3,3) - Jx(1,3)*Jx(3,1))
 198 Jx_inv(3,2) = -(1/detJx) * (Jx(1,1)*Jx(3,2) - Jx(1,2)*Jx(3,1))
 199 Jx_inv(1,3) = +(1/detJx) * (Jx(1,2)*Jx(2,3) - Jx(1,3)*Jx(2,2))
 200 Jx_inv(2,3) = -(1/detJx) * (Jx(1,1)*Jx(2,3) - Jx(1,3)*Jx(2,1))
 201 Jx_inv(3,3) = +(1/detJx) * (Jx(1,1)*Jx(2,2) - Jx(1,2)*Jx(2,1))
 202
 203 Jg = matmul(gradf,Jx_inv)
 204
 205 if(i .eq. 9) then
 206 vol = abs(detJx)*8
 207 u0_tmp = u0*vol
 208 do j = 1,8
 209 tmp = 0
 210 do k = 1,3
 211 tmp = tmp+Jg(j,k)*u0_tmp(k)
 212 end do
 213 tmp_mat(j) = tmp
 214 end do
 215 Floc = Floc-tmp_mat
 216 end if
 217 end do
 218 end if !end for computing Floc
 219
 220 end subroutine hexa
 221
 222 end module module_hexa