1 subroutine elem_search(xl,xts1,xts2,xts3,xsic,elmt, idfile) 2c 3c---------------------------------------------------------------------- 4c This subroutine finds for a particular point in which element it lays 5c and what are its local coordinates 6c 7c input: 8c xl (npro, nshl, nsd) : local node coordinates 9c xts1, xts2, xts3 : point coordinates 10c 11c output: 12c xsic (nsd) : point's local coordinates 13c elmt : element number 14c 15c 16c Elaine Bohr 17c April 2002 18c---------------------------------------------------------------------- 19c 20 use spebc 21 include "common.h" 22 23 dimension shape(nelint,nshl), xsic(nsd), 24 & xl(nelint,nenl,nsd) 25 26 real*8 al(nelint,nenl,nsd), 27 & zi0(nelint,nsd), detaij(nelint), dzi0(nelint,nsd), 28 & neg(nelint), distance(nelint) 29 integer testing(4), found(nelint), negativ(nelint) 30 31 real*8 xts1, xts2, xts3, min 32 integer e, elmt, result, count, counting, find 33 34 tolpt = 0.000001 35 call get_coeff_tet(xl,al) 36 37 detaij(:) = -al(:,2,1)*al(:,3,2)*al(:,4,3) + 38 & al(:,2,1)*al(:,4,2)*al(:,3,3) + al(:,2,2)* 39 & al(:,3,1)*al(:,4,3) - al(:,2,2)*al(:,4,1)* 40 & al(:,3,3) - al(:,2,3)*al(:,3,1)*al(:,4,2)+ 41 & al(:,2,3)*al(:,4,1)*al(:,3,2) 42 43 detaij = 1./detaij 44 45 zi0(:,1) = detaij(:)*((al(:,4,2)*al(:,3,3) 46 & - al(:,3,2)*al(:,4,3))*(xts1-al(:,1,1)) + 47 & (al(:,3,1)*al(:,4,3) 48 & - al(:,4,1)*al(:,3,3))*(xts2-al(:,1,2)) + 49 & (al(:,4,1)*al(:,3,2) 50 & - al(:,3,1)*al(:,4,2))*(xts3-al(:,1,3))) 51 52 53 zi0(:,2) = detaij(:)*((al(:,2,2)*al(:,4,3) 54 & - al(:,4,2)*al(:,2,3))*(xts1-al(:,1,1)) + 55 & (al(:,4,1)*al(:,2,3) 56 & - al(:,2,1)*al(:,4,3))*(xts2-al(:,1,2)) + 57 & (al(:,2,1)*al(:,4,2) 58 & - al(:,4,1)*al(:,2,2))*(xts3-al(:,1,3))) 59 60 zi0(:,3) = detaij(:)*((al(:,3,2)*al(:,2,3) 61 & - al(:,2,2)*al(:,3,3))*(xts1-al(:,1,1)) + 62 & (al(:,2,1)*al(:,3,3) 63 & - al(:,3,1)*al(:,2,3))*(xts2-al(:,1,2)) + 64 & (al(:,3,1)*al(:,2,2) 65 & - al(:,2,1)*al(:,3,2))*(xts3-al(:,1,3))) 66 67c zi0(:,4) = 1 - zi0(:,1) - zi0(:,2) - zi0(:,3) 68 69 elmt = 0 70 counting = 0 71 neg(:) = 0 72 negativ(:) = 0 73 found(:) = 0 74c distance(:) = 0 75 do e = 1, nelint 76 77 count = 0 78 testing(:) = 0 79 if (zi0(e,1).lt.(one+tolpt).and. 80 & zi0(e,1).gt.(zero-tolpt)) then 81 testing(1) = 1 82 count = count + 1 83 endif 84 if (zi0(e,2).lt.(one+tolpt).and. 85 & zi0(e,2).gt.(zero-tolpt)) then 86 testing(2) = 1 87 count = count + 1 88 endif 89 if (zi0(e,3).lt.(one+tolpt).and. 90 & zi0(e,3).gt.(zero-tolpt)) then 91 testing(3) = 1 92 count = count + 1 93 endif 94 if ((1-zi0(e,1)-zi0(e,2)-zi0(e,3)).lt.(one+tolpt).and. 95 & (1-zi0(e,1)-zi0(e,2)-zi0(e,3)).gt.(zero-tolpt)) then 96 testing(4) = 1 97 count = count + 1 98 endif 99 100 result = 1 101 do i = 1, 4 102 result = result*testing(i) 103 enddo 104 105 if (result .eq. 1) then 106 xsic(:) = zi0(e,:) 107 elmt = e 108 109 return 110 111 elseif (count .eq. 3) then 112 counting = counting + 1 113 do i = 1, 3 114 if (testing(i) .eq. 0 .and. 115 & zi0(e,i) .lt. 0.0) then 116 found(counting) = e 117 neg(counting) = zi0(e,i) 118 negativ(counting) = i 119c distance(counting) = sqrt(zi0(e,1)*zi0(e,1)+ 120c & zi0(e,2)*zi0(e,2)+zi0(e,3)*zi0(e,3)) 121 endif 122 enddo 123 if (testing(4) .eq. 0) then 124 found(counting) = e 125 neg(counting) = 1-zi0(e,1)-zi0(e,2)-zi0(e,3) 126 negativ(counting) = 4 127c distance(counting) = sqrt(zi0(e,1)*zi0(e,1)+ 128c & zi0(e,2)*zi0(e,2)+zi0(e,3)*zi0(e,3)) 129 endif 130 endif 131 132 133 enddo 134 135 min = neg(1) 136 elmt = found(1) 137 find = 1 138 do i = 2, counting 139 if (min .lt. neg(i)) then 140 min = neg(i) 141 elmt = found(i) 142 find = i 143 endif 144 enddo 145 if (negativ(find) .eq. 4) then 146 prop = zi0(elmt,1)+zi0(elmt,2)+zi0(elmt,3) 147 xsic(:) = zi0(elmt,:) / prop 148 else 149 xsic(:) = zi0(elmt,:) 150 xsic(negativ(find)) = 0.0 151 endif 152c call error ('elem_search ', 'outrange', idfile) 153 154 return 155 156 157 end 158 159