1c----------------------------------------------------------------------- 2c 3c Natural pressure boundary condition can be calculated with p, the pressure, 4c related (in some prescribed manner) to Q, the flow rate, through the same 5c boundary. To do this efficiently requires us to precompute the integral 6c of N_A over the boundary for each node A and store it in a vector of length 7c nshg (a bit wasteful since only the nodes on the boundary will be non zero 8c in this vector but it is probably slower to index it than to multiply and 9c add the extra zeros....check later). 10c 11c----------------------------------------------------------------------- 12 module pvsQbi 13 14 real*8, allocatable :: NABI(:,:) 15 real*8, allocatable :: NASC(:) 16 integer, allocatable :: ndsurf(:) 17 18 end module 19 20 subroutine finalizeNABI 21 use pvsQbi 22 if( allocated(NABI) ) then 23 deallocate(NABI) 24 endif 25 if( allocated(NASC) ) then 26 deallocate(NASC) 27 endif 28 if( allocated(ndsurf) ) then 29 deallocate(ndsurf) 30 endif 31 end 32 33c----------------------------------------------------------------------- 34c 35c Initialize: 36c 37c----------------------------------------------------------------------- 38 subroutine initNABI( x, shpb ) 39 40 use pointer_data 41 use pvsQbi 42 include "common.h" 43 44 real*8 x(numnp,nsd) 45c 46c use is like 47c 48c NABI=pvsQbi -> NABI 49c 50 dimension shpb(MAXTOP,maxsh,MAXQPT) 51 real*8, allocatable :: tmpshpb(:,:) 52 allocate ( NABI(nshg,3) ) 53 allocate ( NASC(nshg) ) 54 allocate ( ndsurf(nshg) ) 55 56c 57c.... calculate NABI 58c 59 NABI=zero 60 NASC=zero 61 ndsurf=0 62c 63c.... --------------------> boundary elements <-------------------- 64c 65c.... set up parameters 66c 67c intrul = intg (2,itseq) 68c intind = intptb (intrul) 69c 70c.... loop over the boundary elements 71c 72 do iblk = 1, nelblb 73c 74c.... set up the parameters 75c 76 iel = lcblkb(1,iblk) 77 lelCat = lcblkb(2,iblk) 78 lcsyst = lcblkb(3,iblk) 79 iorder = lcblkb(4,iblk) 80 nenl = lcblkb(5,iblk) ! no. of vertices per element 81 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 82 nshl = lcblkb(9,iblk) 83 nshlb = lcblkb(10,iblk) 84 mattyp = lcblkb(7,iblk) 85 ndofl = lcblkb(8,iblk) 86 npro = lcblkb(1,iblk+1) - iel 87 88 89 if(lcsyst.eq.3) lcsyst=nenbl 90c 91 if(lcsyst.eq.3 .or. lcsyst.eq.4) then 92 ngaussb = nintb(lcsyst) 93 else 94 ngaussb = nintb(lcsyst) 95 endif 96 97c 98c.... compute and assemble the residuals corresponding to the 99c boundary integral 100c 101 allocate (tmpshpb(nshl,MAXQPT)) 102 103 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 104 105 call AsBNABI ( x, 106 & tmpshpb, 107 & mienb(iblk)%p, 108 & miBCB(iblk)%p) 109 110 call AsBNASC( x, 111 & tmpshpb, 112 & mienb(iblk)%p, 113 & miBCB(iblk)%p) 114 115 deallocate (tmpshpb) 116 117 enddo 118 119c 120c note that NABI has NOT been communicated. It 121C is the on processor contribution to this vector. It will used to 122C build the on processor contribution to res that will then be made 123C complete via a call to commu. Similarly the LHS usage will create 124C the on-processor contribution to the lhsK. Same for NASC 125c 126 return 127 end 128 129 130