1 subroutine errsmooth(rerr, x, iper, ilwork, 2 & shp, shgl, iBC) 3c 4 use pointer_data 5c 6 include "common.h" 7 include "mpif.h" 8c 9 dimension shp(MAXTOP,maxsh,MAXQPT), 10 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 11 & shpb(MAXTOP,maxsh,MAXQPT), 12 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 13c 14 dimension rerrsm(nshg, 10), rerr(nshg,10), rmass(nshg), 15 & x(nshg,3) 16c 17 dimension ilwork(nlwork), iBC(nshg), iper(nshg) 18 19 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 20 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 21 22c 23c loop over element blocks for the global reconstruction 24c of the smoothed error and lumped mass matrix, rmass 25c 26 rerrsm = zero 27 rmass = zero 28 29 do iblk = 1, nelblk 30c 31c.... set up the parameters 32c 33 nenl = lcblk(5,iblk) ! no. of vertices per element 34 iel = lcblk(1,iblk) 35 lelCat = lcblk(2,iblk) 36 lcsyst = lcblk(3,iblk) 37 iorder = lcblk(4,iblk) 38 nenl = lcblk(5,iblk) ! no. of vertices per element 39 nshl = lcblk(10,iblk) 40 mattyp = lcblk(7,iblk) 41 ndofl = lcblk(8,iblk) 42 nsymdl = lcblk(9,iblk) 43 npro = lcblk(1,iblk+1) - iel 44 ngauss = nint(lcsyst) 45c 46c.... compute and assemble diffusive flux vector residual, qres, 47c and lumped mass matrix, rmass 48 49 allocate (tmpshp(nshl,MAXQPT)) 50 allocate (tmpshgl(nsd,nshl,MAXQPT)) 51 52 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 53 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 54 55 call smooth (rerr, x, 56 & tmpshp, 57 & tmpshgl, 58 & mien(iblk)%p, 59 & rerrsm, 60 & rmass) 61 62 deallocate ( tmpshp ) 63 deallocate ( tmpshgl ) 64 enddo 65c 66 if (numpe > 1) then 67 call commu (rerrsm , ilwork, 10 , 'in ') 68 call commu (rmass , ilwork, 1 , 'in ') 69 endif 70c 71c.... take care of periodic boundary conditions 72c 73 do j= 1,nshg 74 if ((btest(iBC(j),10))) then 75 i = iper(j) 76 rmass(i) = rmass(i) + rmass(j) 77 rerrsm(i,:) = rerrsm(i,:) + rerrsm(j,:) 78 endif 79 enddo 80 81 do j= 1,nshg 82 if ((btest(iBC(j),10))) then 83 i = iper(j) 84 rmass(j) = rmass(i) 85 rerrsm(j,:) = rerrsm(i,:) 86 endif 87 enddo 88c 89c.... invert the diagonal mass matrix and find q 90c 91 rmass = one/rmass 92 93 do i=1, 10 94 rerrsm(:,i) = rmass*rerrsm(:,i) 95 enddo 96 if(numpe > 1) then 97 call commu (rerrsm, ilwork, 10, 'out') 98 endif 99c 100c copy the smoothed error overwriting the original error. 101c 102 103 rerr = rerrsm 104 105 return 106 end 107 108 subroutine smooth (rerr, x, shp, 109 & shgl, ien, 110 & rerrsm, rmass ) 111c 112c---------------------------------------------------------------------- 113c 114c This routine computes and assembles the data corresponding to the 115c interior elements for the global reconstruction of the diffusive 116c flux vector. 117c 118c input: 119c y (nshg,ndof) : Y variables 120c x (numnp,nsd) : nodal coordinates 121c shp (nshape,ngauss) : element shape-functions 122c shgl (nsd,nshape,ngauss) : element local shape-function gradients 123c ien (npro) : nodal connectivity array 124c 125c output: 126c qres (nshg,nflow-1,nsd) : residual vector for diffusive flux 127c rmass (nshg) : lumped mass matrix 128c 129c---------------------------------------------------------------------- 130c 131 include "common.h" 132c 133 dimension rerr(nshg,10), x(numnp,nsd), 134 & shp(nshl,maxsh), 135 & shgl(nsd,nshl,maxsh), 136 & ien(npro,nshl), 137 & rerrsm(nshg,10), rmass(nshg) 138c 139c.... element level declarations 140c 141 dimension rerrl(npro,nshl,10), xl(npro,nenl,nsd), 142 & rerrsml(npro,nshl,10), rmassl(npro,nshl) 143c 144 dimension sgn(npro,nshl), shape(npro,nshl), 145 & shdrv(npro,nsd,nshl), WdetJ(npro), 146 & dxidx(npro,nsd,nsd), shg(npro,nshl,nsd) 147c 148 dimension error(npro,10) 149c 150c.... create the matrix of mode signs for the hierarchic basis 151c functions. 152c 153 if (ipord .gt. 1) then 154 call getsgn(ien,sgn) 155 endif 156c 157c.... gather the variables 158c 159 160 call local(rerr, rerrl, ien, 10, 'gather ') 161 call localx(x, xl, ien, nsd, 'gather ') 162c 163c.... get the element residuals 164c 165 rerrsml = zero 166 rmassl = zero 167 168c 169c.... loop through the integration points 170c 171 172 173 do intp = 1, ngauss 174 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 175c 176c.... create a matrix of shape functions (and derivatives) for each 177c element at this quadrature point. These arrays will contain 178c the correct signs for the hierarchic basis 179c 180 call getshp(shp, shgl, sgn, 181 & shape, shdrv) 182c 183 call e3metric( xl, shdrv, dxidx, 184 & shg, WdetJ) 185 error=zero 186 do n = 1, nshl 187 do i=1,10 188 error(:,i)=error(:,i) + shape(:,n) * rerrl(:,n,i) 189 enddo 190 enddo 191 do i=1,nshl 192 do j=1,10 193 rerrsml(:,i,j) = rerrsml(:,i,j) 194 & + shape(:,i)*WdetJ*error(:,j) 195 enddo 196 197 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 198 enddo 199 200c.... end of the loop over integration points 201c 202 enddo 203c 204c.... assemble the diffusive flux residual 205c 206 call local (rerrsm, rerrsml, ien, 10,'scatter ') 207 call local (rmass, rmassl, ien, 1, 'scatter ') 208c 209 210 return 211 end 212