1 subroutine AsBNABI ( x, shpb, 2 & ienb, iBCB) 3c 4c---------------------------------------------------------------------- 5c 6c This routine computes and assembles the data corresponding to the 7c boundary elements. 8c 9c Zdenek Johan, Winter 1991. (Fortran 90) 10c---------------------------------------------------------------------- 11c 12 use pvsQbi 13 include "common.h" 14c 15 dimension xlb(npro,nenl,nsd), bnorm(npro,nsd), 16 & rl(npro,nshl,nsd), WdetJb(npro), 17 & Wfactor(npro) 18 19 dimension x(numnp,nsd), 20 & shpb(nshl,ngaussb), shglb(nsd,nshl,ngaussb), 21 & ienb(npro,nshl), 22 & iBCB(npro,ndiBCB) 23 24 dimension lnode(27), sgn(npro,nshl), 25 & shpfun(npro,nshl), shdrv(npro,nsd,nshl) 26 27c 28 dimension dxdxib(npro,nsd,nsd), temp(npro), 29 & temp1(npro), temp2(npro), 30 & temp3(npro), 31 & v1(npro,nsd), v2(npro,nsd) 32 33c 34c.... get the matrix of mode signs for the hierarchic basis functions 35c 36 if (ipord .gt. 1) then 37 call getsgn(ienb,sgn) 38 endif 39c 40c.... gather the variables 41c 42 call localx(x, xlb, ienb, nsd, 'gather ') 43c 44c.... get the boundary element residuals 45c 46 rl = zero 47c 48c.... compute the nodes which lie on the boundary (hierarchic) 49c 50 call getbnodes(lnode) 51c 52c.... compute the normal to the boundary. This is achieved by taking 53c the cross product of two vectors in the plane of the 2-d 54c boundary face. 55c 56 if(lcsyst.eq.1) then ! set to curl into element all others out 57 ipt2=2 58 ipt3=3 59 elseif(lcsyst.eq.2) then 60 ipt2=4 61 ipt3=2 62 elseif(lcsyst.eq.3) then 63 ipt2=3 64 ipt3=2 65 elseif(lcsyst.eq.4) then 66 ipt2=2 67 ipt3=4 68 elseif(lcsyst.eq.5) then 69 ipt2=4 70 ipt3=2 71 elseif(lcsyst.eq.6) then 72 ipt2=2 73 ipt3=5 74 endif 75 v1 = xlb(:,ipt2,:) - xlb(:,1,:) 76 v2 = xlb(:,ipt3,:) - xlb(:,1,:) 77c 78c compute cross product 79c 80 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 81 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 82 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 83c 84c mag is area for quads, twice area for tris 85c 86 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 87 bnorm(:,1) = temp1 * temp 88 bnorm(:,2) = temp2 * temp 89 bnorm(:,3) = temp3 * temp 90c 91 if (lcsyst .eq. 1) then 92 Wfactor(:) = one / (four*temp(:)) 93 elseif (lcsyst .eq. 2) then 94 Wfactor(:) = one / (four*temp(:)) 95 elseif (lcsyst .eq. 3) then 96 Wfactor(:) = one / (two*temp(:)) 97 elseif (lcsyst .eq. 4) then 98 Wfactor(:) = one / (four*temp(:)) 99 elseif (lcsyst .eq. 5) then 100 Wfactor(:) = one / (four*temp(:)) 101 elseif (lcsyst .eq. 6) then 102 Wfactor(:) = one / (two*temp(:)) 103 endif 104c 105 do iel=1,npro 106 if (.not.btest(iBCB(iel,1),1)) then 107 Wfactor(iel) = zero ! we want zeros where we are not integrating 108 endif 109 enddo 110c 111 ngaussb = nintb(lcsyst) 112c 113c.... loop through the integration points 114c 115 do intp = 1, ngaussb 116c 117c.... get the hierarchic shape functions at this int point 118c 119 shglb=zero ! protect debugger 120 call getshpb(shpb, shglb, sgn, 121 & shpfun, shdrv) 122c 123 WdetJb(:) = Qwtb(lcsyst,intp) * Wfactor(:) 124c 125c Now lets calculate Integral N_(a:e)^i n_i d Gamma 126c 127 do n = 1, nshlb 128 nodlcl = lnode(n) 129 rl(:,nodlcl,1) = rl(:,nodlcl,1) + 130 & shpfun(:,nodlcl) * bnorm(:,1) * WdetJb(:) 131 rl(:,nodlcl,2) = rl(:,nodlcl,2) + 132 & shpfun(:,nodlcl) * bnorm(:,2) * WdetJb(:) 133 rl(:,nodlcl,3) = rl(:,nodlcl,3) + 134 & shpfun(:,nodlcl) * bnorm(:,3) * WdetJb(:) 135 enddo 136 137 enddo ! quadrature point loop 138c 139c.... assemble the NABI vector 140c 141 call local (NABI, rl, ienb, 3, 'scatter ') 142c 143c push the surf number which we have associated with boundary 144C elements up to the global level in the array ndsurf 145c 146 do iel=1,npro 147 if (iBCB(iel,2) .ne. 0) then 148 iface = iBCB(iel,2) 149 ndsurf(ienb(iel,1:nshlb))=iface 150 endif 151 enddo 152c 153c.... end 154c 155 return 156 end 157 158 159 subroutine AsBNASC ( x, shpb, 160 & ienb, iBCB) 161c 162c---------------------------------------------------------------------- 163c 164c This routine computes and assembles the data corresponding to the 165c boundary elements. 166c 167c Irene Vignon - copied from AsBNABI, Fall 2005. (Fortran 90) 168c---------------------------------------------------------------------- 169c 170 use pvsQbi 171 include "common.h" 172c 173 dimension xlb(npro,nenl,nsd), 174 & rl(npro,nshl), WdetJb(npro), 175 & Wfactor(npro) 176 177 dimension x(numnp,nsd), 178 & shpb(nshl,ngaussb), shglb(nsd,nshl,ngaussb), 179 & ienb(npro,nshl), 180 & iBCB(npro,ndiBCB) 181 182 dimension lnode(27), sgn(npro,nshl), 183 & shpfun(npro,nshl), shdrv(npro,nsd,nshl) 184 185c 186 dimension dxdxib(npro,nsd,nsd), temp(npro), 187 & temp1(npro), temp2(npro), 188 & temp3(npro), 189 & v1(npro,nsd), v2(npro,nsd) 190 191c 192c.... get the matrix of mode signs for the hierarchic basis functions 193c 194 if (ipord .gt. 1) then 195 call getsgn(ienb,sgn) 196 endif 197c 198c.... gather the variables 199c 200 call localx(x, xlb, ienb, nsd, 'gather ') 201c 202c.... get the boundary element residuals 203c 204 rl = zero 205c 206c.... compute the nodes which lie on the boundary (hierarchic) 207c 208 call getbnodes(lnode) 209c 210c.... compute the normal to the boundary. This is achieved by taking 211c the cross product of two vectors in the plane of the 2-d 212c boundary face. 213c 214 if(lcsyst.eq.1) then ! set to curl into element all others out 215 ipt2=2 216 ipt3=3 217 elseif(lcsyst.eq.2) then 218 ipt2=4 219 ipt3=2 220 elseif(lcsyst.eq.3) then 221 ipt2=3 222 ipt3=2 223 elseif(lcsyst.eq.4) then 224 ipt2=2 225 ipt3=4 226 elseif(lcsyst.eq.5) then 227 ipt2=4 228 ipt3=2 229 elseif(lcsyst.eq.6) then 230 ipt2=2 231 ipt3=5 232 endif 233 v1 = xlb(:,ipt2,:) - xlb(:,1,:) 234 v2 = xlb(:,ipt3,:) - xlb(:,1,:) 235c 236c compute cross product 237c 238 temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3) 239 temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3) 240 temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2) 241c 242c mag is area for quads, twice area for tris 243c 244 temp = one / sqrt ( temp1**2 + temp2**2 + temp3**2 ) 245c 246c......here I only want the d Gamma, not the n_i 247 if (lcsyst .eq. 1) then 248 Wfactor(:) = one / (four*temp(:)) 249 elseif (lcsyst .eq. 2) then 250 Wfactor(:) = one / (four*temp(:)) 251 elseif (lcsyst .eq. 3) then 252 Wfactor(:) = one / (two*temp(:)) 253 elseif (lcsyst .eq. 4) then 254 Wfactor(:) = one / (four*temp(:)) 255 elseif (lcsyst .eq. 5) then 256 Wfactor(:) = one / (four*temp(:)) 257 elseif (lcsyst .eq. 6) then 258 Wfactor(:) = one / (two*temp(:)) 259 endif 260c 261 do iel=1,npro 262 if (.not.btest(iBCB(iel,1),1)) then 263 Wfactor(iel) = zero ! we want zeros where we are not integrating 264 endif 265 enddo 266c 267 ngaussb = nintb(lcsyst) 268c 269c.... loop through the integration points 270c 271 do intp = 1, ngaussb 272 273 WdetJb(:) = Qwtb(lcsyst,intp) * Wfactor(:) 274c 275c.... get the hierarchic shape functions at this int point 276c 277 shglb=zero ! protect debugger 278 call getshpb(shpb, shglb, sgn, 279 & shpfun, shdrv) 280 281c 282c Now lets calculate Integral N_(a:e) d Gamma 283c 284 do n = 1, nshlb 285 nodlcl = lnode(n) 286 rl(:,nodlcl) = rl(:,nodlcl) + shpfun(:,nodlcl) * WdetJb(:) 287 enddo 288 289 enddo ! quadrature point loop 290c 291c.... assemble the NASC vector 292c 293 call local (NASC, rl, ienb, 1, 'scatter ') 294c 295c push the surf number which we have associated with boundary 296C elements up to the global level in the array ndsurf 297c 298 do iel=1,npro 299 if (iBCB(iel,2) .ne. 0) then 300 iface = iBCB(iel,2) 301 ndsurf(ienb(iel,1:nshlb))=iface 302 endif 303 enddo 304c 305c.... end 306c 307 return 308 end 309