| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine aspect_ratios(ar)
|
|---|
| 3 | c
|
|---|
| 4 | c 7/6/96
|
|---|
| 5 | c
|
|---|
| 6 | c This routine returns the aspect ratio of a
|
|---|
| 7 | c conglomerate of a set of simplices defined by (tri,nt)
|
|---|
| 8 | c
|
|---|
| 9 | include 'SIZE'
|
|---|
| 10 | include 'GEOM'
|
|---|
| 11 | include 'INPUT'
|
|---|
| 12 | include 'MASS'
|
|---|
| 13 | c
|
|---|
| 14 | real ar(1)
|
|---|
| 15 | real xx(9)
|
|---|
| 16 | c
|
|---|
| 17 | nxyz = lx1*ly1*lz1
|
|---|
| 18 | c
|
|---|
| 19 | if (if3d) then
|
|---|
| 20 | do ie=1,nelv
|
|---|
| 21 | vol = vlsum(bm1(1,1,1,ie),nxyz)
|
|---|
| 22 | x0 = vlsc2(xm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
|
|---|
| 23 | y0 = vlsc2(ym1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
|
|---|
| 24 | z0 = vlsc2(zm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
|
|---|
| 25 | c
|
|---|
| 26 | call rzero(xx,9)
|
|---|
| 27 | do i=1,nxyz
|
|---|
| 28 | x10 = xm1(i,1,1,ie) - x0
|
|---|
| 29 | y10 = ym1(i,1,1,ie) - y0
|
|---|
| 30 | z10 = zm1(i,1,1,ie) - z0
|
|---|
| 31 | c
|
|---|
| 32 | xx(1) = xx(1) + bm1(i,1,1,ie)*x10*x10
|
|---|
| 33 | xx(2) = xx(2) + bm1(i,1,1,ie)*x10*y10
|
|---|
| 34 | xx(3) = xx(3) + bm1(i,1,1,ie)*x10*z10
|
|---|
| 35 | xx(4) = xx(4) + bm1(i,1,1,ie)*y10*x10
|
|---|
| 36 | xx(5) = xx(5) + bm1(i,1,1,ie)*y10*y10
|
|---|
| 37 | xx(6) = xx(6) + bm1(i,1,1,ie)*y10*z10
|
|---|
| 38 | xx(7) = xx(7) + bm1(i,1,1,ie)*z10*x10
|
|---|
| 39 | xx(8) = xx(8) + bm1(i,1,1,ie)*z10*y10
|
|---|
| 40 | xx(9) = xx(9) + bm1(i,1,1,ie)*z10*z10
|
|---|
| 41 | enddo
|
|---|
| 42 | vi = 1./vol
|
|---|
| 43 | call cmult(xx,vi,9)
|
|---|
| 44 | c call eig3(xx,eign,eig1)
|
|---|
| 45 | call eig2(xx,eign,eig1)
|
|---|
| 46 | ar(ie) = sqrt(eign/eig1)
|
|---|
| 47 | enddo
|
|---|
| 48 | else
|
|---|
| 49 | do ie=1,nelv
|
|---|
| 50 | vol = vlsum(bm1(1,1,1,ie),nxyz)
|
|---|
| 51 | x0 = vlsc2(xm1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
|
|---|
| 52 | y0 = vlsc2(ym1(1,1,1,ie),bm1(1,1,1,ie),nxyz)/vol
|
|---|
| 53 | c
|
|---|
| 54 | call rzero(xx,4)
|
|---|
| 55 | do i=1,nxyz
|
|---|
| 56 | x10 = xm1(i,1,1,ie) - x0
|
|---|
| 57 | y10 = ym1(i,1,1,ie) - y0
|
|---|
| 58 | xx(1) = xx(1) + bm1(i,1,1,ie)*x10*x10
|
|---|
| 59 | xx(2) = xx(2) + bm1(i,1,1,ie)*x10*y10
|
|---|
| 60 | xx(3) = xx(3) + bm1(i,1,1,ie)*y10*x10
|
|---|
| 61 | xx(4) = xx(4) + bm1(i,1,1,ie)*y10*y10
|
|---|
| 62 | enddo
|
|---|
| 63 | vi = 1./vol
|
|---|
| 64 | call cmult(xx,vi,4)
|
|---|
| 65 | call eig2(xx,eign,eig1)
|
|---|
| 66 | c write(6,6) ie,vol,eign,eig1
|
|---|
| 67 | c 6 format(i5,' veig:',1p3e16.6)
|
|---|
| 68 | ar(ie) = sqrt(eign/eig1)
|
|---|
| 69 | enddo
|
|---|
| 70 | endif
|
|---|
| 71 | c
|
|---|
| 72 | c do ie=1,nelv
|
|---|
| 73 | c write(6,*) ar(ie),ie,' aspect'
|
|---|
| 74 | c enddo
|
|---|
| 75 | c
|
|---|
| 76 | return
|
|---|
| 77 | end
|
|---|
| 78 | c-----------------------------------------------------------------------
|
|---|
| 79 | subroutine eig2(AA,eign,eig1)
|
|---|
| 80 | c
|
|---|
| 81 | c return max and min eigenvalues of a 2x2 matrix
|
|---|
| 82 | c
|
|---|
| 83 | real aa(2,2)
|
|---|
| 84 | c
|
|---|
| 85 | a = aa(1,1)
|
|---|
| 86 | b = aa(1,2)
|
|---|
| 87 | c = aa(2,1)
|
|---|
| 88 | d = aa(2,2)
|
|---|
| 89 | c
|
|---|
| 90 | qa = 1.
|
|---|
| 91 | qb = -(a+d)
|
|---|
| 92 | qc = (a*d - c*b)
|
|---|
| 93 | call quadratic_h(eign,eig1,qa,qb,qc)
|
|---|
| 94 | c
|
|---|
| 95 | return
|
|---|
| 96 | end
|
|---|
| 97 | c-----------------------------------------------------------------------
|
|---|
| 98 | subroutine quadratic_h(x1,x2,a,b,c)
|
|---|
| 99 | c
|
|---|
| 100 | c Stable routine for computation of real roots of quadratic
|
|---|
| 101 | c
|
|---|
| 102 | c The "_h" denotes the hack below so we don't need to worry
|
|---|
| 103 | c about complex arithmetic in the near-double root case. pff 1/22/97
|
|---|
| 104 | c
|
|---|
| 105 | c Upon return, | x1 | >= | x2 |
|
|---|
| 106 | c
|
|---|
| 107 | x1 = 0.
|
|---|
| 108 | x2 = 0.
|
|---|
| 109 | c
|
|---|
| 110 | if (a.eq.0.) then
|
|---|
| 111 | if (b.eq.0) then
|
|---|
| 112 | write(6,10) x1,x2,a,b,c
|
|---|
| 113 | return
|
|---|
| 114 | endif
|
|---|
| 115 | x1 = -c/b
|
|---|
| 116 | write(6,11) x1,a,b,c
|
|---|
| 117 | return
|
|---|
| 118 | endif
|
|---|
| 119 | c
|
|---|
| 120 | d = b*b - 4.*a*c
|
|---|
| 121 | if (d.lt.0) then
|
|---|
| 122 | c write(6,12) a,b,c,d
|
|---|
| 123 | c hack, for this application we'll assume d<0 by epsilon, and just set
|
|---|
| 124 | d = 0
|
|---|
| 125 | endif
|
|---|
| 126 | if (d.gt.0) d = sqrt(d)
|
|---|
| 127 | c
|
|---|
| 128 | q = -0.5 * (b + b/abs(b) * d)
|
|---|
| 129 | x1 = q/a
|
|---|
| 130 | x2 = c/q
|
|---|
| 131 | c
|
|---|
| 132 | if (abs(x2).gt.abs(x1)) then
|
|---|
| 133 | tmp = x2
|
|---|
| 134 | x2 = x1
|
|---|
| 135 | x1 = tmp
|
|---|
| 136 | endif
|
|---|
| 137 | c
|
|---|
| 138 | 10 format('ERROR: Both a & b zero in routine quadratic NO ROOTS.'
|
|---|
| 139 | $ ,1p5e12.4)
|
|---|
| 140 | 11 format('ERROR: a = 0 in routine quadratic, only one root.'
|
|---|
| 141 | $ ,1p5e12.4)
|
|---|
| 142 | 12 format('ERROR: negative discriminate in routine quadratic.'
|
|---|
| 143 | $ ,1p5e12.4)
|
|---|
| 144 | c
|
|---|
| 145 | return
|
|---|
| 146 | end
|
|---|
| 147 | c-----------------------------------------------------------------------
|
|---|
| 148 | subroutine out_sem(iel)
|
|---|
| 149 | c
|
|---|
| 150 | include 'SIZE'
|
|---|
| 151 | include 'INPUT'
|
|---|
| 152 | include 'GEOM'
|
|---|
| 153 | c
|
|---|
| 154 | c
|
|---|
| 155 | open(unit=33,file='v1')
|
|---|
| 156 | if (iel.eq.0) then
|
|---|
| 157 | do ie=1,nelv
|
|---|
| 158 | write(33,33) xm1( 1, 1,1,ie),ym1( 1, 1,1,ie)
|
|---|
| 159 | write(33,33) xm1(lx1, 1,1,ie),ym1(lx1, 1,1,ie)
|
|---|
| 160 | write(33,33) xm1(lx1,ly1,1,ie),ym1(lx1,ly1,1,ie)
|
|---|
| 161 | write(33,33) xm1( 1,ly1,1,ie),ym1( 1,ly1,1,ie)
|
|---|
| 162 | enddo
|
|---|
| 163 | else
|
|---|
| 164 | ie = iel
|
|---|
| 165 | write(33,33) xm1( 1, 1,1,ie),ym1( 1, 1,1,ie)
|
|---|
| 166 | write(33,33) xm1(lx1, 1,1,ie),ym1(lx1, 1,1,ie)
|
|---|
| 167 | write(33,33) xm1(lx1,ly1,1,ie),ym1(lx1,ly1,1,ie)
|
|---|
| 168 | write(33,33) xm1( 1,ly1,1,ie),ym1( 1,ly1,1,ie)
|
|---|
| 169 | endif
|
|---|
| 170 | 33 format(f14.6)
|
|---|
| 171 | close(unit=33)
|
|---|
| 172 | return
|
|---|
| 173 | end
|
|---|
| 174 | c-----------------------------------------------------------------------
|
|---|
| 175 | subroutine gradm11(ux,uy,uz,u,e)
|
|---|
| 176 | c
|
|---|
| 177 | c Compute gradient of T -- mesh 1 to mesh 1 (vel. to vel.)
|
|---|
| 178 | c
|
|---|
| 179 | c Single element case
|
|---|
| 180 | c
|
|---|
| 181 | include 'SIZE'
|
|---|
| 182 | include 'DXYZ'
|
|---|
| 183 | include 'GEOM'
|
|---|
| 184 | include 'INPUT'
|
|---|
| 185 | include 'TSTEP'
|
|---|
| 186 | c
|
|---|
| 187 | parameter (lxyz=lx1*ly1*lz1)
|
|---|
| 188 | real ux(lxyz),uy(lxyz),uz(lxyz),u(lxyz,1)
|
|---|
| 189 | c
|
|---|
| 190 | common /ctmp1/ ur(lxyz),us(lxyz),ut(lxyz)
|
|---|
| 191 | c
|
|---|
| 192 | integer e
|
|---|
| 193 |
|
|---|
| 194 | c
|
|---|
| 195 | N = lx1-1
|
|---|
| 196 | if (if3d) then
|
|---|
| 197 | call local_grad3(ur,us,ut,u,N,e,dxm1,dxtm1)
|
|---|
| 198 | do i=1,lxyz
|
|---|
| 199 | ux(i) = jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
|
|---|
| 200 | $ + us(i)*sxm1(i,1,1,e)
|
|---|
| 201 | $ + ut(i)*txm1(i,1,1,e) )
|
|---|
| 202 | uy(i) = jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
|
|---|
| 203 | $ + us(i)*sym1(i,1,1,e)
|
|---|
| 204 | $ + ut(i)*tym1(i,1,1,e) )
|
|---|
| 205 | uz(i) = jacmi(i,e)*(ur(i)*rzm1(i,1,1,e)
|
|---|
| 206 | $ + us(i)*szm1(i,1,1,e)
|
|---|
| 207 | $ + ut(i)*tzm1(i,1,1,e) )
|
|---|
| 208 | enddo
|
|---|
| 209 | else
|
|---|
| 210 | if (ifaxis) call setaxdy (ifrzer(e))
|
|---|
| 211 | call local_grad2(ur,us,u,N,e,dxm1,dytm1)
|
|---|
| 212 | do i=1,lxyz
|
|---|
| 213 | ux(i) =jacmi(i,e)*(ur(i)*rxm1(i,1,1,e)
|
|---|
| 214 | $ + us(i)*sxm1(i,1,1,e) )
|
|---|
| 215 | uy(i) =jacmi(i,e)*(ur(i)*rym1(i,1,1,e)
|
|---|
| 216 | $ + us(i)*sym1(i,1,1,e) )
|
|---|
| 217 | enddo
|
|---|
| 218 | endif
|
|---|
| 219 | c
|
|---|
| 220 | return
|
|---|
| 221 | end
|
|---|
| 222 | c-----------------------------------------------------------------------
|
|---|
