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