femwind/fortran/module_ndt_assembly.f90

   1 module module_ndt_assembly
   2        use module_hexa
   3
   4 contains
   5
   6 subroutine ndt_assembly(                              &
   7     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
   8     ifms, ifme, kfms,kfme, jfms, jfme,            &
   9     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  10     ifts, ifte, kfts, kfte, jfts,jfte,            &
  11     A, X,Y,Z, iflags,                          & !Input from femwind, U0, V0, W0 not used in hexa to construct Kloc or JG
  12     K)                                             !Global stiffness matrix output
  13
  14 implicit none
  15
  16
  17 !*** arguments
  18
  19 integer, intent(in)::                     &
  20     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions in elements
  21     ifms, ifme, kfms,kfme, jfms, jfme,            &
  22     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  23     ifts, ifte, kfts, kfte, jfts,jfte
  24
  25
  26
  27
  28 integer, parameter:: msize=14
  29 real, intent(in), dimension(3,3):: A
  30 real, intent(in), dimension(ifms:ifme, kfms:kfme, jfms: jfme):: X,Y,Z!spatial grid
  31 !Input for hexa
  32 integer, intent(in)::iflags
  33
  34 real, intent(out), dimension(ifms:ifme, kfms:kfme, jfms:jfme,1:msize)::K
  35
  36
  37 !*** local
  38
  39 integer:: ie1, ie2, ie3, ic1, ic2, ic3, iloc, i
  40 real ::  Kloc(8,8), Jg(8,3), Xloc(3,8), vol
  41
  42 !*** u0loc is an input for module_hexa, but is not used to construct K. Do I need to define this?
  43 !*** integer, dimension(3,1,1), save ::iflags = reshape((/1,0,1/),(/3,1,1/)) !define iflags to construct JG and Kloc in hexa
  44
  45
  46
  47 !** executable
  48 Xloc = 99999.
  49 K = 0.
  50
  51 do ie2=jfts,jfte
  52     do ie3=kfts, kfte
  53         do ie1=ifts, ifte
  54             do ic2=0,1
  55                 do ic3=0,1
  56                     do ic1=0,1
  57                         iloc=1+ic1+2*(ic2+2*ic3);  !local index of the node in the element
  58                         Xloc(1,iloc)=X(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  59                         Xloc(2,iloc)=Y(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  60                         Xloc(3,iloc)=Z(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  61                     enddo
  62                 enddo
  63             enddo
  64             !write(*,*)'Xloc='
  65             !write(*,'((8f9.4))')((Xloc(i,iloc),iloc=1,8),i=1,3)
  66
  67             call hexa(A,Xloc,Kloc,Jg,vol,iflags)
  68
  69             !write(*,*)'Kloc='
  70             !write(*,'((8f9.4))')((Kloc(i,iloc),iloc=1,8),i=1,8)
  71
  72          !print *, 'Kloc(1,1)=',Kloc(1,1)
  73          !print *, ie1, ie2, ie3, K(1,1,2,1)
  74             K(ie1  ,ie3  ,ie2  , 1) =   K(ie1  ,ie3  ,ie2  , 1) + Kloc( 1,  1)
  75             K(ie1  ,ie3  ,ie2  , 2) =   K(ie1  ,ie3  ,ie2  , 2) + Kloc( 1,  2)
  76             K(ie1  ,ie3  ,ie2  , 4) =   K(ie1  ,ie3  ,ie2  , 4) + Kloc( 1,  3)
  77             K(ie1  ,ie3  ,ie2  , 5) =   K(ie1  ,ie3  ,ie2  , 5) + Kloc( 1,  4)
  78             K(ie1  ,ie3  ,ie2  ,10) =   K(ie1  ,ie3  ,ie2  ,10) + Kloc( 1,  5)
  79             K(ie1  ,ie3  ,ie2  ,11) =   K(ie1  ,ie3  ,ie2  ,11) + Kloc( 1,  6)
  80             K(ie1  ,ie3  ,ie2  ,13) =   K(ie1  ,ie3  ,ie2  ,13) + Kloc( 1,  7)
  81             K(ie1  ,ie3  ,ie2  ,14) =   K(ie1  ,ie3  ,ie2  ,14) + Kloc( 1,  8)
  82             K(ie1+1,ie3  ,ie2  , 1) =   K(ie1+1,ie3  ,ie2  , 1) + Kloc( 2,  2)
  83             K(ie1+1,ie3  ,ie2  , 3) =   K(ie1+1,ie3  ,ie2  , 3) + Kloc( 2,  3)
  84             K(ie1+1,ie3  ,ie2  , 4) =   K(ie1+1,ie3  ,ie2  , 4) + Kloc( 2,  4)
  85             K(ie1+1,ie3  ,ie2  , 9) =   K(ie1+1,ie3  ,ie2  , 9) + Kloc( 2,  5)
  86             K(ie1+1,ie3  ,ie2  ,10) =   K(ie1+1,ie3  ,ie2  ,10) + Kloc( 2,  6)
  87             K(ie1+1,ie3  ,ie2  ,12) =   K(ie1+1,ie3  ,ie2  ,12) + Kloc( 2,  7)
  88             K(ie1+1,ie3  ,ie2  ,13) =   K(ie1+1,ie3  ,ie2  ,13) + Kloc( 2,  8)
  89             K(ie1  ,ie3  ,ie2+1, 1) =   K(ie1  ,ie3  ,ie2+1, 1) + Kloc( 3,  3)
  90             K(ie1  ,ie3  ,ie2+1, 2) =   K(ie1  ,ie3  ,ie2+1, 2) + Kloc( 3,  4)
  91             K(ie1  ,ie3  ,ie2+1, 7) =   K(ie1  ,ie3  ,ie2+1, 7) + Kloc( 3,  5)
  92             K(ie1  ,ie3  ,ie2+1, 8) =   K(ie1  ,ie3  ,ie2+1, 8) + Kloc( 3,  6)
  93             K(ie1  ,ie3  ,ie2+1,10) =   K(ie1  ,ie3  ,ie2+1,10) + Kloc( 3,  7)
  94             K(ie1  ,ie3  ,ie2+1,11) =   K(ie1  ,ie3  ,ie2+1,11) + Kloc( 3,  8)
  95             K(ie1+1,ie3  ,ie2+1, 1) =   K(ie1+1,ie3  ,ie2+1, 1) + Kloc( 4,  4)
  96             K(ie1+1,ie3  ,ie2+1, 6) =   K(ie1+1,ie3  ,ie2+1, 6) + Kloc( 4,  5)
  97             K(ie1+1,ie3  ,ie2+1, 7) =   K(ie1+1,ie3  ,ie2+1, 7) + Kloc( 4,  6)
  98             K(ie1+1,ie3  ,ie2+1, 9) =   K(ie1+1,ie3  ,ie2+1, 9) + Kloc( 4,  7)
  99             K(ie1+1,ie3  ,ie2+1,10) =   K(ie1+1,ie3  ,ie2+1,10) + Kloc( 4,  8)
 100             K(ie1  ,ie3+1,ie2  , 1) =   K(ie1  ,ie3+1,ie2  , 1) + Kloc( 5,  5)
 101             K(ie1  ,ie3+1,ie2  , 2) =   K(ie1  ,ie3+1,ie2  , 2) + Kloc( 5,  6)
 102             K(ie1  ,ie3+1,ie2  , 4) =   K(ie1  ,ie3+1,ie2  , 4) + Kloc( 5,  7)
 103             K(ie1  ,ie3+1,ie2  , 5) =   K(ie1  ,ie3+1,ie2  , 5) + Kloc( 5,  8)
 104             K(ie1+1,ie3+1,ie2  , 1) =   K(ie1+1,ie3+1,ie2  , 1) + Kloc( 6,  6)
 105             K(ie1+1,ie3+1,ie2  , 3) =   K(ie1+1,ie3+1,ie2  , 3) + Kloc( 6,  7)
 106             K(ie1+1,ie3+1,ie2  , 4) =   K(ie1+1,ie3+1,ie2  , 4) + Kloc( 6,  8)
 107             K(ie1  ,ie3+1,ie2+1, 1) =   K(ie1  ,ie3+1,ie2+1, 1) + Kloc( 7,  7)
 108             K(ie1  ,ie3+1,ie2+1, 2) =   K(ie1  ,ie3+1,ie2+1, 2) + Kloc( 7,  8)
 109             K(ie1+1,ie3+1,ie2+1, 1) =   K(ie1+1,ie3+1,ie2+1, 1) + Kloc( 8,  8)
 110        !print *, ie1, ie2, ie3, K(1,1,2,1)
 111         enddo
 112     enddo
 113 enddo
 114
 115 end subroutine ndt_assembly
 116 end module  module_ndt_assembly
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136