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 20c----------------------------------------------------------------------- 21c 22c Initialize: 23c 24c----------------------------------------------------------------------- 25 subroutine initNABI( x, shpb ) 26 27 use pointer_data 28 use pvsQbi 29 include "common.h" 30 31 real*8 x(numnp,nsd) 32c 33c use is like 34c 35c NABI=pvsQbi -> NABI 36c 37 dimension shpb(MAXTOP,maxsh,MAXQPT) 38 real*8, allocatable :: tmpshpb(:,:) 39 allocate ( NABI(nshg,3) ) 40 allocate ( NASC(nshg) ) 41 allocate ( ndsurf(nshg) ) 42 43c 44c.... calculate NABI 45c 46 NABI=zero 47 NASC=zero 48 ndsurf=0 49c 50c.... --------------------> boundary elements <-------------------- 51c 52c.... set up parameters 53c 54c intrul = intg (2,itseq) 55c intind = intptb (intrul) 56c 57c.... loop over the boundary elements 58c 59 do iblk = 1, nelblb 60c 61c.... set up the parameters 62c 63 iel = lcblkb(1,iblk) 64 lelCat = lcblkb(2,iblk) 65 lcsyst = lcblkb(3,iblk) 66 iorder = lcblkb(4,iblk) 67 nenl = lcblkb(5,iblk) ! no. of vertices per element 68 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 69 nshl = lcblkb(9,iblk) 70 nshlb = lcblkb(10,iblk) 71 mattyp = lcblkb(7,iblk) 72 ndofl = lcblkb(8,iblk) 73 npro = lcblkb(1,iblk+1) - iel 74 75 76 if(lcsyst.eq.3) lcsyst=nenbl 77c 78 if(lcsyst.eq.3 .or. lcsyst.eq.4) then 79 ngaussb = nintb(lcsyst) 80 else 81 ngaussb = nintb(lcsyst) 82 endif 83 84c 85c.... compute and assemble the residuals corresponding to the 86c boundary integral 87c 88 allocate (tmpshpb(nshl,MAXQPT)) 89 90 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 91 92 call AsBNABI ( x, 93 & tmpshpb, 94 & mienb(iblk)%p, 95 & miBCB(iblk)%p) 96 97 call AsBNASC( x, 98 & tmpshpb, 99 & mienb(iblk)%p, 100 & miBCB(iblk)%p) 101 102 deallocate (tmpshpb) 103 104 enddo 105 106c 107c note that NABI has NOT been communicated. It 108C is the on processor contribution to this vector. It will used to 109C build the on processor contribution to res that will then be made 110C complete via a call to commu. Similarly the LHS usage will create 111C the on-processor contribution to the lhsK. Same for NASC 112c 113 return 114 end 115 116 117