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), Floc(8), Jg(8,3)
  41 real ::  Xloc(3,8), u0loc(3)
  42
  43 !*** u0loc is an input for module_hexa, but is not used to construct K. Do I need to define this?
  44 !*** integer, dimension(3,1,1), save ::iflags = reshape((/1,0,1/),(/3,1,1/)) !define iflags to construct JG and Kloc in hexa
  45
  46
  47
  48 !** executable
  49 Xloc = 99999.
  50 K = 0.
  51
  52 do ie2=jfts,jfte
  53     do ie3=kfts, kfte
  54         do ie1=ifts, ifte
  55             do ic2=0,1
  56                 do ic3=0,1
  57                     do ic1=0,1
  58                         iloc=1+ic1+2*(ic2+2*ic3);  !local index of the node in the element
  59                         Xloc(1,iloc)=X(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  60                         Xloc(2,iloc)=Y(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  61                         Xloc(3,iloc)=Z(ie1 + ic1, ie3 + ic3, ie2 + ic2)
  62                     enddo
  63                 enddo
  64             enddo
  65             !write(*,*)'Xloc='
  66             !write(*,'((8f9.4))')((Xloc(i,iloc),iloc=1,8),i=1,3)
  67
  68             call hexa(A,Xloc,u0loc,Kloc,Floc,Jg,iflags)
  69
  70             !write(*,*)'Kloc='
  71             !write(*,'((8f9.4))')((Kloc(i,iloc),iloc=1,8),i=1,8)
  72
  73          !print *, 'Kloc(1,1)=',Kloc(1,1)
  74          !print *, ie1, ie2, ie3, K(1,1,2,1)
  75             K(ie1  ,ie3  ,ie2  , 1) =   K(ie1  ,ie3  ,ie2  , 1) + Kloc( 1,  1)
  76             K(ie1  ,ie3  ,ie2  , 2) =   K(ie1  ,ie3  ,ie2  , 2) + Kloc( 1,  2)
  77             K(ie1  ,ie3  ,ie2  , 4) =   K(ie1  ,ie3  ,ie2  , 4) + Kloc( 1,  3)
  78             K(ie1  ,ie3  ,ie2  , 5) =   K(ie1  ,ie3  ,ie2  , 5) + Kloc( 1,  4)
  79             K(ie1  ,ie3  ,ie2  ,10) =   K(ie1  ,ie3  ,ie2  ,10) + Kloc( 1,  5)
  80             K(ie1  ,ie3  ,ie2  ,11) =   K(ie1  ,ie3  ,ie2  ,11) + Kloc( 1,  6)
  81             K(ie1  ,ie3  ,ie2  ,13) =   K(ie1  ,ie3  ,ie2  ,13) + Kloc( 1,  7)
  82             K(ie1  ,ie3  ,ie2  ,14) =   K(ie1  ,ie3  ,ie2  ,14) + Kloc( 1,  8)
  83             K(ie1+1,ie3  ,ie2  , 1) =   K(ie1+1,ie3  ,ie2  , 1) + Kloc( 2,  2)
  84             K(ie1+1,ie3  ,ie2  , 3) =   K(ie1+1,ie3  ,ie2  , 3) + Kloc( 2,  3)
  85             K(ie1+1,ie3  ,ie2  , 4) =   K(ie1+1,ie3  ,ie2  , 4) + Kloc( 2,  4)
  86             K(ie1+1,ie3  ,ie2  , 9) =   K(ie1+1,ie3  ,ie2  , 9) + Kloc( 2,  5)
  87             K(ie1+1,ie3  ,ie2  ,10) =   K(ie1+1,ie3  ,ie2  ,10) + Kloc( 2,  6)
  88             K(ie1+1,ie3  ,ie2  ,12) =   K(ie1+1,ie3  ,ie2  ,12) + Kloc( 2,  7)
  89             K(ie1+1,ie3  ,ie2  ,13) =   K(ie1+1,ie3  ,ie2  ,13) + Kloc( 2,  8)
  90             K(ie1  ,ie3  ,ie2+1, 1) =   K(ie1  ,ie3  ,ie2+1, 1) + Kloc( 3,  3)
  91             K(ie1  ,ie3  ,ie2+1, 2) =   K(ie1  ,ie3  ,ie2+1, 2) + Kloc( 3,  4)
  92             K(ie1  ,ie3  ,ie2+1, 7) =   K(ie1  ,ie3  ,ie2+1, 7) + Kloc( 3,  5)
  93             K(ie1  ,ie3  ,ie2+1, 8) =   K(ie1  ,ie3  ,ie2+1, 8) + Kloc( 3,  6)
  94             K(ie1  ,ie3  ,ie2+1,10) =   K(ie1  ,ie3  ,ie2+1,10) + Kloc( 3,  7)
  95             K(ie1  ,ie3  ,ie2+1,11) =   K(ie1  ,ie3  ,ie2+1,11) + Kloc( 3,  8)
  96             K(ie1+1,ie3  ,ie2+1, 1) =   K(ie1+1,ie3  ,ie2+1, 1) + Kloc( 4,  4)
  97             K(ie1+1,ie3  ,ie2+1, 6) =   K(ie1+1,ie3  ,ie2+1, 6) + Kloc( 4,  5)
  98             K(ie1+1,ie3  ,ie2+1, 7) =   K(ie1+1,ie3  ,ie2+1, 7) + Kloc( 4,  6)
  99             K(ie1+1,ie3  ,ie2+1, 9) =   K(ie1+1,ie3  ,ie2+1, 9) + Kloc( 4,  7)
 100             K(ie1+1,ie3  ,ie2+1,10) =   K(ie1+1,ie3  ,ie2+1,10) + Kloc( 4,  8)
 101             K(ie1  ,ie3+1,ie2  , 1) =   K(ie1  ,ie3+1,ie2  , 1) + Kloc( 5,  5)
 102             K(ie1  ,ie3+1,ie2  , 2) =   K(ie1  ,ie3+1,ie2  , 2) + Kloc( 5,  6)
 103             K(ie1  ,ie3+1,ie2  , 4) =   K(ie1  ,ie3+1,ie2  , 4) + Kloc( 5,  7)
 104             K(ie1  ,ie3+1,ie2  , 5) =   K(ie1  ,ie3+1,ie2  , 5) + Kloc( 5,  8)
 105             K(ie1+1,ie3+1,ie2  , 1) =   K(ie1+1,ie3+1,ie2  , 1) + Kloc( 6,  6)
 106             K(ie1+1,ie3+1,ie2  , 3) =   K(ie1+1,ie3+1,ie2  , 3) + Kloc( 6,  7)
 107             K(ie1+1,ie3+1,ie2  , 4) =   K(ie1+1,ie3+1,ie2  , 4) + Kloc( 6,  8)
 108             K(ie1  ,ie3+1,ie2+1, 1) =   K(ie1  ,ie3+1,ie2+1, 1) + Kloc( 7,  7)
 109             K(ie1  ,ie3+1,ie2+1, 2) =   K(ie1  ,ie3+1,ie2+1, 2) + Kloc( 7,  8)
 110             K(ie1+1,ie3+1,ie2+1, 1) =   K(ie1+1,ie3+1,ie2+1, 1) + Kloc( 8,  8)
 111        !print *, ie1, ie2, ie3, K(1,1,2,1)
 112         enddo
 113     enddo
 114 enddo
 115
 116 end subroutine ndt_assembly
 117 end module  module_ndt_assembly
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137