1 2 subroutine timeseries(ycl, xl, ien, sgn) 3 4 use timedataC 5 include "common.h" 6 7 dimension shape(nshl), ycl(npro,nshl,ndofl), 8 & ien(npro,nshl), xl(npro,nenl,nsd), 9 & sgn(npro,nshl) 10 real*8 al(npro,nenl,nsd), 11 & zi0(npro,nsd), detaij(npro), dzi0(npro,nsd), 12 & m11(npro), m12(npro), m13(npro), m21(npro), m22(npro), 13 & m23(npro), m31(npro), m32(npro), m33(npro), 14 & r1(npro), r2(npro), r3(npro), shgradl(nshl,nsd) 15 16 real*8 xts1, xts2, xts3 17 real*8 soln(ndof) 18 integer e, founde 19 20 do jj = 1, ntspts 21 founde = 0 22 if(statptts(jj,1).gt.0) then 23 if(statptts(jj,1).eq.iblkts) then 24 if(lcsyst.eq.2) then ! hex 25 call shphex (ipord, parptts(jj,:),shape(:), 26 & shgradl(:,:)) 27 elseif(lcsyst.eq.1) then 28 call shptet (ipord, parptts(jj,:),shape(:), 29 & shgradl(:,:)) 30 endif 31 founde=statptts(jj,2) 32 endif 33 else 34 xts1 = ptts(jj,1) 35 xts2 = ptts(jj,2) 36 xts3 = ptts(jj,3) 37 38 if(lcsyst.eq.2) then ! hex 39 40 call get_a_not_hex(xl,al) ! get mapping poly. coeff. 41 42c... get initial guess for Newton Iteration procedure 43 44 detaij(:) = -al(:,2,1)*al(:,3,2)*al(:,4,3) + 45 & al(:,2,1)*al(:,4,2)*al(:,3,3) + al(:,2,2)* 46 & al(:,3,1)*al(:,4,3) - al(:,2,2)*al(:,4,1)* 47 & al(:,3,3) - al(:,2,3)*al(:,3,1)*al(:,4,2)+ 48 & al(:,2,3)*al(:,4,1)*al(:,3,2) 49 50 detaij = 1./detaij 51 52 zi0(:,1) = detaij(:)*((al(:,4,2)*al(:,3,3) 53 & - al(:,3,2)*al(:,4,3))*(xts1-al(:,1,1)) + 54 & (al(:,3,1)*al(:,4,3) 55 & - al(:,4,1)*al(:,3,3))*(xts2-al(:,1,2)) + 56 & (al(:,4,1)*al(:,3,2) 57 & - al(:,3,1)*al(:,4,2))*(xts3-al(:,1,3))) 58 59 60 zi0(:,2) = detaij(:)*((al(:,2,2)*al(:,4,3) 61 & - al(:,4,2)*al(:,2,3))*(xts1-al(:,1,1)) + 62 & (al(:,4,1)*al(:,2,3) 63 & - al(:,2,1)*al(:,4,3))*(xts2-al(:,1,2)) + 64 & (al(:,2,1)*al(:,4,2) 65 & - al(:,4,1)*al(:,2,2))*(xts3-al(:,1,3))) 66 67 zi0(:,3) = detaij(:)*((al(:,3,2)*al(:,2,3) 68 & - al(:,2,2)*al(:,3,3))*(xts1-al(:,1,1)) + 69 & (al(:,2,1)*al(:,3,3) 70 & - al(:,3,1)*al(:,2,3))*(xts2-al(:,1,2)) + 71 & (al(:,3,1)*al(:,2,2) 72 & - al(:,2,1)*al(:,3,2))*(xts3-al(:,1,3))) 73 74 75c... iterate to convergence 76 77 do it = 1, iterat 78 79c... build matrix 80 81 m11(:)=al(:,2,1)+al(:,5,1)*zi0(:,2)+al(:,7,1)*zi0(:,3) 82 & +al(:,8,1)*zi0(:,2)*zi0(:,3) 83 m12(:)=al(:,3,1)+al(:,5,1)*zi0(:,1)+al(:,6,1)*zi0(:,3) 84 & +al(:,8,1)*zi0(:,1)*zi0(:,3) 85 m13(:)=al(:,4,1)+al(:,6,1)*zi0(:,2)+al(:,7,1)*zi0(:,1) 86 & +al(:,8,1)*zi0(:,1)*zi0(:,2) 87 88 m21(:)=al(:,2,2)+al(:,5,2)*zi0(:,2)+al(:,7,2)*zi0(:,3) 89 & +al(:,8,2)*zi0(:,2)*zi0(:,3) 90 m22(:)=al(:,3,2)+al(:,5,2)*zi0(:,1)+al(:,6,2)*zi0(:,3) 91 & +al(:,8,2)*zi0(:,1)*zi0(:,3) 92 m23(:)=al(:,4,2)+al(:,6,2)*zi0(:,2)+al(:,7,2)*zi0(:,1) 93 & +al(:,8,2)*zi0(:,1)*zi0(:,2) 94 95 m31(:)=al(:,2,3)+al(:,5,3)*zi0(:,2)+al(:,7,3)*zi0(:,3) 96 & +al(:,8,3)*zi0(:,2)*zi0(:,3) 97 m32(:)=al(:,3,3)+al(:,5,3)*zi0(:,1)+al(:,6,3)*zi0(:,3) 98 & +al(:,8,3)*zi0(:,1)*zi0(:,3) 99 m33(:)=al(:,4,3)+al(:,6,3)*zi0(:,2)+al(:,7,3)*zi0(:,1) 100 & +al(:,8,3)*zi0(:,1)*zi0(:,2) 101 102 103c... build rhs 104 105 r1(:)=al(:,1,1)+al(:,2,1)*zi0(:,1)+al(:,3,1)*zi0(:,2)+ 106 & al(:,4,1)*zi0(:,3)+al(:,5,1)*zi0(:,1)*zi0(:,2)+ 107 & al(:,6,1)*zi0(:,2)*zi0(:,3)+al(:,7,1)* 108 & zi0(:,1)*zi0(:,3)+al(:,8,1)*zi0(:,1)* 109 & zi0(:,2)*zi0(:,3) - xts1 110 111 r2(:)=al(:,1,2)+al(:,2,2)*zi0(:,1)+al(:,3,2)*zi0(:,2)+ 112 & al(:,4,2)*zi0(:,3)+al(:,5,2)*zi0(:,1)*zi0(:,2)+ 113 & al(:,6,2)*zi0(:,2)*zi0(:,3)+al(:,7,2)* 114 & zi0(:,1)*zi0(:,3)+al(:,8,2)*zi0(:,1)* 115 & zi0(:,2)*zi0(:,3) - xts2 116 117 r3(:)=al(:,1,3)+al(:,2,3)*zi0(:,1)+al(:,3,3)*zi0(:,2)+ 118 & al(:,4,3)*zi0(:,3)+al(:,5,3)*zi0(:,1)*zi0(:,2)+ 119 & al(:,6,3)*zi0(:,2)*zi0(:,3)+al(:,7,3)* 120 & zi0(:,1)*zi0(:,3)+al(:,8,3)*zi0(:,1)* 121 & zi0(:,2)*zi0(:,3) - xts3 122 123c... get solution 124 125 detaij = m11*m22*m33-m11*m23*m32-m21*m12*m33+ 126 & m21*m13*m32+m31*m12*m23-m31*m13*m22 127 128 detaij = 1./detaij 129 130 dzi0(:,1) = -detaij(:)*((m22(:)*m33(:)-m23(:)*m32(:))* 131 & r1(:) + (m13(:)*m32(:)-m12(:)*m33(:))*r2(:) + 132 & (m12(:)*m23(:)-m13(:)*m22(:))*r3(:)) 133 dzi0(:,2) = -detaij(:)*((m23(:)*m31(:)-m21(:)*m32(:))* 134 & r1(:) + (m11(:)*m33(:)-m13(:)*m31(:))*r2(:) + 135 & (m13(:)*m21(:)-m11(:)*m23(:))*r3(:)) 136 dzi0(:,3) = -detaij(:)*((m21(:)*m32(:)-m22(:)*m31(:))* 137 & r1(:) + (m12(:)*m31(:)-m11(:)*m32(:))*r2(:) + 138 & (m11(:)*m22(:)-m12(:)*m21(:))*r3(:)) 139 140 zi0(:,:) = zi0(:,:) + dzi0(:,:) 141 142 enddo 143 144 do e = 1, npro 145 if ((abs(zi0(e,1)).lt.(one+tolpt)).and. 146 & (abs(zi0(e,2)).lt.(one+tolpt)).and. 147 & (abs(zi0(e,3)).lt.(one+tolpt))) then ! got the element 148 149 call shphex (ipord, zi0(e,:),shape(:), 150 & shgradl(:,:)) 151 152 founde=e 153 exit 154 endif 155 156 enddo 157 elseif (lcsyst.eq.1) then !tet 158 159 call get_a_not_tet(xl,al) 160 161c 162c solve for r, s, t for each elements 163c 164 do e = 1, npro 165 detaij(e) = al(e,2,1)*(-al(e,3,2)*al(e,4,3) + 166 & al(e,4,2)*al(e,3,3)) + al(e,2,2)* 167 & (al(e,3,1)*al(e,4,3) - al(e,4,1)* 168 & al(e,3,3)) + al(e,2,3)*(-al(e,3,1)*al(e,4,2)+ 169 & al(e,4,1)*al(e,3,2)) 170 171 detaij(e) = 1./detaij(e) 172 173 zi0(e,1) = detaij(e)*((al(e,4,2)*al(e,3,3) 174 & - al(e,3,2)*al(e,4,3))*(xts1-al(e,1,1)) + 175 & (al(e,3,1)*al(e,4,3) 176 & - al(e,4,1)*al(e,3,3))*(xts2-al(e,1,2)) + 177 & (al(e,4,1)*al(e,3,2) 178 & - al(e,3,1)*al(e,4,2))*(xts3-al(e,1,3))) 179 180 zi0(e,2) = detaij(e)*((al(e,2,2)*al(e,4,3) 181 & - al(e,4,2)*al(e,2,3))*(xts1-al(e,1,1)) + 182 & (al(e,4,1)*al(e,2,3) 183 & - al(e,2,1)*al(e,4,3))*(xts2-al(e,1,2)) + 184 & (al(e,2,1)*al(e,4,2) 185 & - al(e,4,1)*al(e,2,2))*(xts3-al(e,1,3))) 186 187 zi0(e,3) = detaij(e)*((al(e,3,2)*al(e,2,3) 188 & - al(e,2,2)*al(e,3,3))*(xts1-al(e,1,1)) + 189 & (al(e,2,1)*al(e,3,3) 190 & - al(e,3,1)*al(e,2,3))*(xts2-al(e,1,2)) + 191 & (al(e,3,1)*al(e,2,2) 192 & - al(e,2,1)*al(e,3,2))*(xts3-al(e,1,3))) 193 194 if ((zi0(e,1)+zi0(e,2)+zi0(e,3)).lt.(one+tolpt).and. 195 & zi0(e,1).lt.(one +tolpt).and. !should not be necessary; the limit in the tet is already defined by the previous line 196 & zi0(e,1).gt.(zero-tolpt).and. 197 & zi0(e,2).lt.(one +tolpt).and. !should not be necessary 198 & zi0(e,2).gt.(zero-tolpt).and. 199 & zi0(e,3).lt.(one +tolpt).and. !should not be necessary 200 & zi0(e,3).gt.(zero-tolpt)) then 201 202 call shptet (ipord, zi0(e,:), shape(:), 203 & shgradl(:,:)) 204 205 founde=e 206 exit 207 endif 208 enddo 209 endif 210 211 if(founde.ne.0) then !store values for next time steps 212 statptts(jj,1)=iblkts 213 statptts(jj,2)=founde 214 parptts(jj,:)=zi0(founde,:) 215 else 216 if(iblkts.eq.nelblk) then !not found in last elm-blk of this part 217 statptts(jj,1)=nelblk+1 ! not searched for following steps 218 endif 219 endif 220 endif 221 222 if(founde.ne.0) then 223 soln(1:ndofl) = zero 224 225 do i = 1,nenl 226 soln(1:ndofl) = soln(1:ndofl) 227 & +ycl(founde,i,1:ndofl)*shape(i) 228 enddo 229 do i = 1+nenl,nshl 230 soln(1:ndofl) = soln(1:ndofl) 231 & +ycl(founde,i,1:ndofl)*shape(i)*sgn(founde,i) 232 enddo 233 234 varts(jj,:) = soln(:) 235 endif 236 237 enddo 238 239 return 240 end 241 242