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