femwind/fortran/ndt_assembly_test.f90

   1 program ndt_assembly_test
   2
   3 use module_ndt_assembly
   4 use module_hexa
   5 use module_utils
   6
   7 implicit none
   8
   9 real:: Amat(3,3)
  10 real, pointer:: Xmat(:,:,:),Ymat(:,:,:), Zmat(:,:,:), Kmat(:,:,:,:)
  11 real, pointer::a1(:), a2(:), a3(:), a4(:), X(:,:,:), Y(:,:,:), Z(:,:,:), Km(:,:,:,:)
  12 integer ::n1(2),n2(3),n3(3), n4(3), m(3)
  13
  14 integer :: msize, &
  15     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  16     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  17     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  18     ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  19     iats, iate, jats, jate, iams,iame,                        &         ! Amat bounds
  20     jams, jame
  21
  22 integer :: i,j,k,jx
  23 integer :: iflags = 1                !Flags to construct K in hexa module
  24 integer:: ksize(4)                               ! Global Stifness Matrix Dimensions in Matlab
  25
  26 call read_array_nd(a1,n1,'A') !Recovering X-Matrix and dimension of X matrix
  27 if (n1(1).ne.3.or.n1(2).ne.3)then
  28     call crash('A must be 3 by 3')
  29     stop
  30 endif
  31
  32 Amat = reshape(a1,n1)
  33
  34 ! read input arrays in ikj index ordering and tight bounds
  35 call read_array(X,'X')        !Recovering X-Matrix and dimension of X matrix
  36 call read_array(Y,'Y')        !Recovering X-Matrix and dimension of X matrix
  37 call read_array(Z,'Z')        !Recovering X-Matrix and dimension of X matrix
  38
  39 n2 = shape(X)
  40
  41 ifts = 1
  42 ifte = n2(1)-1
  43 jfts = 1
  44 jfte = n2(3)-1
  45 kfts = 1
  46 kfte = n2(2)-1
  47 msize = 14
  48
  49 ifms = ifts
  50 ifme = ifte+2
  51 jfms = jfts
  52 jfme = jfte+2
  53 kfms = kfts
  54 kfme = kfte+2
  55
  56
  57 allocate(Kmat(ifms:ifme,kfms:kfme,jfms:jfme, 1:msize))
  58
  59 allocate(Xmat(ifms:ifme,kfms:kfme,jfms:jfme))
  60 allocate(Ymat(ifms:ifme,kfms:kfme,jfms:jfme))
  61 allocate(Zmat(ifms:ifme,kfms:kfme,jfms:jfme))
  62
  63 ! copy the input data to tile sized bounds
  64 do j=jfts,jfte+1
  65   do k=kfts,kfte+1
  66     do i=ifts,ifte+1
  67         Xmat(i,k,j) = X(i,k,j)
  68         Ymat(i,k,j) = Y(i,k,j)
  69         Zmat(i,k,j) = Z(i,k,j)
  70     enddo
  71   enddo
  72 enddo
  73
  74 !write(*,'(a)')'calling ndt_assembly'
  75 call ndt_assembly(  &
  76   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  77   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  78   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  79   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  80   Amat,Xmat,Ymat,Zmat, iflags, Kmat)
  81
  82 !write(*,'(a,3i8)')'copying the output data to array size ',n2,msize
  83 allocate(Km(ifts:ifte+1,kfts:kfte+1,jfts:jfte+1, 1:msize))
  84 do j=jfts,jfte+1
  85   do k=kfts,kfte+1
  86     do i=ifts,ifte+1
  87         do jx = 1,msize
  88             Km(i,k,j,jx)=Kmat(i,k,j,jx)
  89         enddo
  90     enddo
  91   enddo
  92 enddo
  93
  94
  95 ksize = (/ifte-ifts+2,kfte-kfts+2,jfte-jfts+2,msize/)
  96 call write_array_nd(reshape(Km,(/product(ksize)/)),ksize,'K')
  97 !***print *,(/product(ksize)/)
  98
  99
 100 end program ndt_assembly_test