1 subroutine solvecon(y, x, iBC, BC, 2 & iper, ilwork, shp, shgl) 3 4c--------------------------------------------------------------------- 5c This subroutine is to calculate the constarint for the redistancing 6c scalar of the level set method. This is to prevent interface from 7c moving by applying the condition that the volume must stay constant 8c in each element when the redisatnce step is applied. 9c-------------------------------------------------------------------- 10c 11c 12 use pointer_data 13c 14 include "common.h" 15 include "mpif.h" 16 include "auxmpi.h" 17c 18 dimension y(nshg,ndof), 19 & x(numnp,nsd), iBC(nshg), 20 & BC(nshg,ndofBC), ilwork(nlwork), 21 & iper(nshg) 22c 23c 24 dimension shp(MAXTOP,maxsh,MAXQPT), 25 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 26 & v_lambda(nshg), hprime(nshg), 27 & v_lambda1(nshg), v_lambda2(nshg), 28 & rmass(nshg) 29c 30 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 31 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 32 33c 34c ... intialize 35c 36 rmass = zero 37 v_lambda = zero 38 v_lambda1 = zero 39 v_lambda2 = zero 40 hprime = zero 41c 42c ... loop over element blocks 43c 44 do iblk = 1, nelblk 45c 46c.... set up the parameters 47c 48 nenl = lcblk(5,iblk) ! no. of vertices per element 49 iel = lcblk(1,iblk) 50 lelCat = lcblk(2,iblk) 51 lcsyst = lcblk(3,iblk) 52 iorder = lcblk(4,iblk) 53 nenl = lcblk(5,iblk) ! no. of vertices per element 54 nshl = lcblk(10,iblk) 55 mattyp = lcblk(7,iblk) 56 ndofl = lcblk(8,iblk) 57 nsymdl = lcblk(9,iblk) 58 npro = lcblk(1,iblk+1) - iel 59 ngauss = nint(lcsyst) 60c 61c.... compute and assemble the constarint factor, and mass matrix 62c 63 64 allocate (tmpshp(nshl,MAXQPT)) 65 allocate (tmpshgl(nsd,nshl,MAXQPT)) 66 67 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 68 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 69c 70 call volcon (y, x, tmpshp, 71 & tmpshgl, mien(iblk)%p, rmass, 72 & v_lambda1, hprime, v_lambda2) 73 74 deallocate ( tmpshp ) 75 deallocate ( tmpshgl ) 76 enddo 77 78c 79c ... multiple processor communication 80c 81 if (numpe > 1) then 82 call commu (v_lambda1 , ilwork, 1 , 'in ') 83 call commu (v_lambda2 , ilwork, 1 , 'in ') 84 call commu (hprime , ilwork, 1 , 'in ') 85 call commu (rmass , ilwork, 1 , 'in ') 86 endif 87c 88c.... take care of periodic boundary conditions 89c 90 do j= 1,nshg 91 if (btest(iBC(j),10)) then 92 i = iper(j) 93 rmass(i) = rmass(i) + rmass(j) 94 v_lambda1(i) = v_lambda1(i) + v_lambda1(j) 95 v_lambda2(i) = v_lambda2(i) + v_lambda2(j) 96 hprime(i) = hprime(i) + hprime(j) 97 endif 98 enddo 99c 100 do j= 1,nshg 101 if (btest(iBC(j),10)) then 102 i = iper(j) 103 rmass(j) = rmass(i) 104 v_lambda1(j) = v_lambda1(i) 105 v_lambda2(j) = v_lambda2(i) 106 hprime(j) = hprime(i) 107 endif 108 enddo 109c 110c ... calculation of constraint factor 111c 112 rmass = one/rmass 113 v_lambda1 = v_lambda1*rmass ! numerator of lambda 114 v_lambda2 = v_lambda2*rmass ! denominator of lambda 115 v_lambda = v_lambda1/(v_lambda2+epsM**2) 116 hprime = hprime*rmass 117 v_lambda = v_lambda*hprime 118c 119c ... commu out for the multiple processor 120c 121 if(numpe > 1) then 122 call commu (v_lambda, ilwork, 1, 'out') 123 endif 124c 125c ... the following commented lines are for the different way of getting 126c the denominator of constraint (lambda) calculation 127c$$$ hprime=zero 128c$$$ do kk=1, nshg 129c$$$ if (abs (y(kk,6)) .le. epsilon_ls) then 130c$$$ hprime(kk) = (0.5/epsilon_ls) * (1 131c$$$ & + cos(pi*y(kk,6)/epsilon_ls)) 132c$$$ endif 133c$$$ enddo 134c$$$ y(:,7) = y(:,7)+v_lambda*hprime/dtgl 135c 136c ... the vlome constraint applied on the second scalar 137c 138 y(:,7) = y(:,7)+v_lambda/dtgl 139c 140 return 141 end 142c 143c 144c 145 146 subroutine volcon (y, x, shp, 147 & shgl, ien, rmass, 148 & v_lambda1, hprime, v_lambda2) 149 150c--------------------------------------------------------------------- 151c 152c This subroutine is to calculate the element contribution to the 153c constraint factor and mass matrix. 154c 155c--------------------------------------------------------------------- 156 include "common.h" 157c 158 dimension y(nshg,ndof), x(numnp,nsd), 159 & shp(nshl,maxsh), 160 & shgl(nsd,nshl,maxsh), 161 & ien(npro,nshl), 162 & qres(nshg,idflx), rmass(nshg) 163c 164c.... element level declarations 165c 166 dimension ycl(npro,nshl,ndof), xl(npro,nenl,nsd), 167 & rmassl(npro,nshl) 168 dimension sgn(npro,nshape), v_lambdal1(npro,nshl), 169 & v_lambda1(nshg), hprimel(npro,nshl), 170 & hprime(nshg), v_lambdal2(npro,nshl), 171 & v_lambda2(nshg) 172c 173c local arrays 174c 175 dimension shg(npro,nshl,nsd), 176 & dxidx(npro,nsd,nsd), WdetJ(npro) 177c 178 dimension shape(npro,nshl), 179 & shdrv(npro,nsd,nshl) 180c 181c 182c.... for volume constraint calculation of redistancing step 183c 184 dimension Sclr(npro), Sclrtmp(npro), 185 & h_prime(npro), tmp1(npro), 186 & tmp2(npro), 187 & v_lambdatmp(npro), v_lambdal(npro,nshl) 188c$$$ & ,hprimel(npro,nshl), v_lambdal1(npro,nshl), 189c$$$ & v_lambdal2(npro,nshl) 190c above arrays must be uncommented for alternate method included below (commented) 191 real epsilon_tmp 192 193c 194c.... create the matrix of mode signs for the hierarchic basis 195c functions. 196c 197 if (ipord .gt. 1) then 198 call getsgn(ien,sgn) 199 endif 200c 201c.... gather the variables 202c 203 204 call localy(y, ycl, ien, ndof, 'gather ') 205 call localx(x, xl, ien, nsd, 'gather ') 206c 207c.... get the element contributions of the numerator and denominator 208c of lambda 209c 210 rmassl = zero 211 v_lambdal1= zero 212 v_lambdal2= zero 213 hprimel = zero 214c 215c.... loop through the integration points 216c 217 do intp = 1, ngauss 218 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 219c 220c.... create a matrix of shape functions (and derivatives) for each 221c element at this quadrature point. These arrays will contain 222c the correct signs for the hierarchic basis 223c 224 call getshp(shp, shgl, sgn, 225 & shape, shdrv) 226c 227c.... initialize 228c 229 h_prime = zero 230 sclr = zero 231 sclrtmp = zero 232 v_lambdatmp=zero 233 tmp1 =zero 234 tmp2 =zero 235 236c 237c.... ---------------------> Element Metrics <----------------------- 238c 239 call e3metric( xl, shdrv, dxidx, 240 & shg, WdetJ) 241 242c 243 do i = 1, nshl 244c 245c y(intp)=SUM_{a=1}^nshl (N_a(intp) Ya) 246c 247 Sclr = Sclr + shape(:,i) * ycl(:,i,7) !d^kbar 248 sclrtmp = sclrtmp + shape(:,i) * ycl(:,i,6) !d^0 249 enddo 250 251 if (isclr .eq. 2) then 252 epsilon_tmp = epsilon_lsd 253 else 254 epsilon_tmp = epsilon_ls 255 endif 256 257 do i=1,npro 258 if (abs (Sclrtmp(i)) .le. epsilon_tmp) then 259 h_prime(i) = (0.5/epsilon_tmp) * (1 260 & + cos(pi*Sclrtmp(i)/epsilon_tmp)) 261 tmp1(i)=-h_prime(i)*(sclr(i)-sclrtmp(i))*dtgl 262 tmp2(i)=h_prime(i)**2 263c v_lambdatmp(i)=tmp1(i)/(tmp2(i)+ epsM) 264 endif 265 enddo 266c 267 do i=1,nshl 268 v_lambdal1(:,i)=v_lambdal1(:,i)+ shape(:,i)*WdetJ*tmp1 269 v_lambdal2(:,i)=v_lambdal2(:,i)+ shape(:,i)*WdetJ*tmp2 270c v_lambdal(:,i)=v_lambdal(:,i)+ shape(:,i)*WdetJ*v_lambdatmp 271 hprimel(:,i) = hprimel(:,i) + shape(:,i)*WdetJ*h_prime 272 rmassl(:,i) =rmassl(:,i) + shape(:,i)*WdetJ 273 enddo 274c 275c.... end of the loop over integration points 276c 277 enddo 278c 279 call local (v_lambda1, v_lambdal1, ien, 1, 'scatter ') 280 call local (v_lambda2, v_lambdal2, ien, 1, 'scatter ') 281 call local (hprime, hprimel, ien, 1, 'scatter ') 282 call local (rmass, rmassl, ien, 1, 'scatter ') 283c 284 return 285 end 286 287