| c------------------------------------------------------------------ |
| subroutine setuppc |
| c------------------------------------------------------------------ |
| c Generate diagonal preconditioner for CG. |
| c Preconditioner computed in this subroutine is correct only |
| c for collocation point in element interior, on conforming face |
| c interior and conforming edge. |
| c------------------------------------------------------------------ |
| |
| include 'header.h' |
| |
| double precision dxtm1_2(lx1,lx1), rdtime |
| integer ie,k,i,j,q,isize |
| |
| do j=1,lx1 |
| do i=1,lx1 |
| dxtm1_2(i,j)=dxtm1(i,j)**2 |
| end do |
| end do |
| |
| rdtime=1.d0/dtime |
| |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ie,isize,i,j,k,q) |
| do ie = 1, nelt |
| call r_init(dpcelm(1,1,1,ie),nxyz,0.d0) |
| isize=size_e(ie) |
| do k = 1, lx1 |
| do j = 1, lx1 |
| do i = 1, lx1 |
| do q = 1, lx1 |
| dpcelm(i,j,k,ie) = dpcelm(i,j,k,ie) + |
| & g1m1_s(q,j,k,isize) * dxtm1_2(i,q) + |
| & g1m1_s(i,q,k,isize) * dxtm1_2(j,q) + |
| & g1m1_s(i,j,q,isize) * dxtm1_2(k,q) |
| end do |
| dpcelm(i,j,k,ie)=visc*dpcelm(i,j,k,ie)+ |
| & rdtime*bm1_s(i,j,k,isize) |
| end do |
| end do |
| end do |
| end do |
| c$OMP END PARALLEL DO |
| |
| c.....do the stiffness summation |
| call dssum |
| |
| c.....take inverse. |
| |
| call reciprocal(dpcelm,ntot) |
| |
| c.....compute preconditioner on mortar points. NOTE: dpcmor for |
| c nonconforming cases will be corrected in subroutine setpcmo |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I) |
| do i=1,nmor |
| dpcmor(i)=1.d0/dpcmor(i) |
| end do |
| c$OMP END PARALLEL DO |
| |
| return |
| end |
| |
| |
| c-------------------------------------------------------------- |
| subroutine setpcmo_pre |
| c-------------------------------------------------------------- |
| c pre-compute elemental contribution to preconditioner |
| c for all situations |
| c-------------------------------------------------------------- |
| |
| include 'header.h' |
| |
| integer element_size, i, j, ii, jj, col |
| double precision |
| & p(lx1,lx1,lx1), p0(lx1,lx1,lx1), mtemp(lx1,lx1), |
| & temp(lx1,lx1,lx1), temp1(lx1,lx1), tmp(lx1,lx1),tig(lx1) |
| |
| c.....corners on face of type 3 |
| |
| call r_init(tcpre,lx1*lx1,0.d0) |
| call r_init(tmp,lx1*lx1,0.d0) |
| call r_init(tig,5,0.d0) |
| tig(1) =1.d0 |
| tmp(1,1) =1.d0 |
| |
| c.....tcpre results from mapping a unit spike field (unity at |
| c collocation point (1,1), zero elsewhere) on an entire element |
| c face to the (1,1) segment of a nonconforming face |
| do i=2,lx1-1 |
| do j=1,lx1 |
| tmp(i,1) = tmp(i,1)+ qbnew(i-1,j,1)*tig(j) |
| end do |
| end do |
| |
| do col=1,lx1 |
| tcpre(col,1)=tmp(col,1) |
| |
| do j=2,lx1-1 |
| do i=1,lx1 |
| tcpre(col,j) = tcpre(col,j) + qbnew(j-1,i,1)* |
| & tmp(col,i) |
| end do |
| end do |
| end do |
| |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(element_size,i,j,p,temp, |
| c$OMP& mtemp,temp1,p0,ii,jj) |
| do element_size=1,refine_max |
| |
| c.......for conforming cases |
| |
| c.......pcmor_c (i,j,element_size) records the intermediate value |
| c (preconditioner=1/pcmor_c) of the preconditor on collocation |
| c point (i,j) on a conforming face of an element of size |
| c element_size. |
| |
| do j=1,lx1/2+1 |
| do i=j,lx1/2+1 |
| call r_init(p,nxyz,0.d0) |
| p(i,j,1)=1.d0 |
| call laplacian(temp,p,element_size) |
| pcmor_c(i,j,element_size)=temp(i,j,1) |
| pcmor_c(lx1+1-i,j,element_size)=temp(i,j,1) |
| pcmor_c(j,i,element_size)=temp(i,j,1) |
| pcmor_c(lx1+1-j,i,element_size)=temp(i,j,1) |
| pcmor_c(j,lx1+1-i,element_size)=temp(i,j,1) |
| pcmor_c(lx1+1-j,lx1+1-i,element_size)=temp(i,j,1) |
| pcmor_c(i,lx1+1-j,element_size)=temp(i,j,1) |
| pcmor_c(lx1+1-i,lx1+1-j,element_size)=temp(i,j,1) |
| end do |
| end do |
| |
| c.......for nonconforming cases |
| |
| c.......nonconforming face interior |
| |
| c.......pcmor_nc1(i,j,ii,jj,element_size) records the intermediate |
| c preconditioner value on collocation point (i,j) on mortar |
| c (ii,jj) on a nonconforming face of an element of size element_ |
| c size |
| do j=2,lx1 |
| do i=j,lx1 |
| call r_init(mtemp,lx1*lx1,0.d0) |
| call r_init(p,nxyz,0.d0) |
| mtemp(i,j)=1.d0 |
| c...........when i, j=lx1, mortar points are duplicated, so mtemp needs |
| c to be doubled. |
| if(i.eq.lx1)mtemp(i,j)=mtemp(i,j)*2.d0 |
| if(j.eq.lx1)mtemp(i,j)=mtemp(i,j)*2.d0 |
| call transf_nc(mtemp,p) |
| call laplacian(temp,p,element_size) |
| call transfb_nc1(temp1,temp) |
| |
| c...........values at points (i,j) and (j,i) are the same |
| pcmor_nc1(i,j,1,1,element_size)=temp1(i,j) |
| pcmor_nc1(j,i,1,1,element_size)=temp1(i,j) |
| end do |
| |
| c.........when i, j=lx1, mortar points are duplicated. so pcmor_nc1 needs |
| c to be doubled on those points |
| pcmor_nc1(lx1,j,1,1,element_size)= |
| & pcmor_nc1(lx1,j,1,1,element_size)*2.d0 |
| pcmor_nc1(j,lx1,1,1,element_size)= |
| & pcmor_nc1(lx1,j,1,1,element_size) |
| |
| end do |
| pcmor_nc1(lx1,lx1,1,1,element_size)= |
| & pcmor_nc1(lx1,lx1,1,1,element_size)*2.d0 |
| |
| c.......nonconforming edges |
| j=1 |
| do i=2,lx1 |
| call r_init(mtemp,lx1*lx1,0.d0) |
| call r_init(p,nxyz,0.d0) |
| call r_init(p0,nxyz,0.d0) |
| mtemp(i,j)=1.d0 |
| if(i.eq.lx1)mtemp(i,j)=2.d0 |
| call transf_nc(mtemp,p) |
| call laplacian(temp,p,element_size) |
| call transfb_nc1(temp1,temp) |
| pcmor_nc1(i,j,1,1,element_size)=temp1(i,j) |
| pcmor_nc1(j,i,1,1,element_size)=temp1(i,j) |
| do ii=1,lx1 |
| c...........p0 is for the case that a nonconforming edge is shared by |
| c two conforming faces |
| p0(ii,1,1)=p(ii,1,1) |
| do jj=1,lx1 |
| c.............now p is for the case that a nonconforming edge is shared |
| c by nonconforming faces |
| p(ii,1,jj)=p(ii,jj,1) |
| end do |
| end do |
| |
| call laplacian(temp,p,element_size) |
| call transfb_nc2(temp1,temp) |
| |
| c.........pcmor_nc2(i,j,ii,jj,element_size) gives the intermediate |
| c preconditioner value on collocation point (i,j) on a |
| c nonconforming face of an element with size size_element |
| |
| pcmor_nc2(i,j,1,1,element_size)=temp1(i,j)*2.d0 |
| pcmor_nc2(j,i,1,1,element_size)= |
| & pcmor_nc2(i,j,1,1,element_size) |
| |
| call laplacian(temp,p0,element_size) |
| call transfb_nc0(temp1,temp) |
| |
| c.........pcmor_nc0(i,j,ii,jj,element_size) gives the intermediate |
| c preconditioner value on collocation point (i,j) on a |
| c conforming face of an element, which shares a nonconforming |
| c edge with another conforming face |
| pcmor_nc0(i,j,1,1,element_size)=temp1(i,j) |
| pcmor_nc0(j,i,1,1,element_size)=temp1(i,j) |
| end do |
| pcmor_nc1(lx1,j,1,1,element_size)= |
| & pcmor_nc1(lx1,j,1,1,element_size)*2.d0 |
| pcmor_nc1(j,lx1,1,1,element_size)= |
| & pcmor_nc1(lx1,j,1,1,element_size) |
| pcmor_nc2(lx1,j,1,1,element_size)= |
| & pcmor_nc2(lx1,j,1,1,element_size)*2.d0 |
| pcmor_nc2(j,lx1,1,1,element_size)= |
| & pcmor_nc2(lx1,j,1,1,element_size) |
| pcmor_nc0(lx1,j,1,1,element_size)= |
| & pcmor_nc0(lx1,j,1,1,element_size)*2.d0 |
| pcmor_nc0(j,lx1,1,1,element_size)= |
| & pcmor_nc0(lx1,j,1,1,element_size) |
| |
| c.......symmetrical copy |
| do i=1,lx1-1 |
| pcmor_nc1(i,j,1,2,element_size)= |
| & pcmor_nc1(lx1+1-i,j,1,1,element_size) |
| pcmor_nc0(i,j,1,2,element_size)= |
| & pcmor_nc0(lx1+1-i,j,1,1,element_size) |
| pcmor_nc2(i,j,1,2,element_size)= |
| & pcmor_nc2(lx1+1-i,j,1,1,element_size) |
| end do |
| |
| do j=2,lx1 |
| do i=1,lx1-1 |
| pcmor_nc1(i,j,1,2,element_size)= |
| & pcmor_nc1(lx1+1-i,j,1,1,element_size) |
| end do |
| i=lx1 |
| pcmor_nc1(i,j,1,2,element_size)= |
| & pcmor_nc1(lx1+1-i,j,1,1,element_size) |
| pcmor_nc0(i,j,1,2,element_size)= |
| & pcmor_nc0(lx1+1-i,j,1,1,element_size) |
| pcmor_nc2(i,j,1,2,element_size)= |
| & pcmor_nc2(lx1+1-i,j,1,1,element_size) |
| end do |
| |
| j=1 |
| i=1 |
| pcmor_nc1(i,j,2,1,element_size)= |
| & pcmor_nc1(i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,1,element_size)= |
| & pcmor_nc0(i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,1,element_size)= |
| & pcmor_nc2(i,lx1+1-j,1,1,element_size) |
| do j=2,lx1-1 |
| i=1 |
| pcmor_nc1(i,j,2,1,element_size)= |
| & pcmor_nc1(i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,1,element_size)= |
| & pcmor_nc0(i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,1,element_size)= |
| & pcmor_nc2(i,lx1+1-j,1,1,element_size) |
| do i=2,lx1 |
| pcmor_nc1(i,j,2,1,element_size)= |
| & pcmor_nc1(i,lx1+1-j,1,1,element_size) |
| end do |
| end do |
| |
| j=lx1 |
| do i=2,lx1 |
| pcmor_nc1(i,j,2,1,element_size)= |
| & pcmor_nc1(i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,1,element_size)= |
| & pcmor_nc0(i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,1,element_size)= |
| & pcmor_nc2(i,lx1+1-j,1,1,element_size) |
| end do |
| |
| j=1 |
| i=lx1 |
| pcmor_nc1(i,j,2,2,element_size)= |
| & pcmor_nc1(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,2,element_size)= |
| & pcmor_nc0(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,2,element_size)= |
| & pcmor_nc2(lx1+1-i,lx1+1-j,1,1,element_size) |
| |
| do j=2,lx1-1 |
| do i=2,lx1-1 |
| pcmor_nc1(i,j,2,2,element_size)= |
| & pcmor_nc1(lx1+1-i,lx1+1-j,1,1,element_size) |
| end do |
| i=lx1 |
| pcmor_nc1(i,j,2,2,element_size)= |
| & pcmor_nc1(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,2,element_size)= |
| & pcmor_nc0(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,2,element_size)= |
| & pcmor_nc2(lx1+1-i,lx1+1-j,1,1,element_size) |
| end do |
| j=lx1 |
| do i=2,lx1-1 |
| pcmor_nc1(i,j,2,2,element_size)= |
| & pcmor_nc1(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc0(i,j,2,2,element_size)= |
| & pcmor_nc0(lx1+1-i,lx1+1-j,1,1,element_size) |
| pcmor_nc2(i,j,2,2,element_size)= |
| & pcmor_nc2(lx1+1-i,lx1+1-j,1,1,element_size) |
| end do |
| |
| |
| c.......vertices shared by at least one nonconforming face or edge |
| |
| c.......Among three edges and three faces sharing a vertex on an element |
| c situation 1: only one edge is nonconforming |
| c situation 2: two edges are nonconforming |
| c situation 3: three edges are nonconforming |
| c situation 4: one face is nonconforming |
| c situation 5: one face and one edge are nonconforming |
| c situation 6: two faces are nonconforming |
| c situation 7: three faces are nonconforming |
| |
| call r_init(p0,nxyz,0.d0) |
| p0(1,1,1)=1.d0 |
| call laplacian(temp,p0,element_size) |
| pcmor_cor(8,element_size)=temp(1,1,1) |
| |
| c.......situation 1 |
| call r_init(p0,nxyz,0.d0) |
| do i=1,lx1 |
| p0(i,1,1)=tcpre(i,1) |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_e(1,pcmor_cor(1,element_size),temp) |
| |
| c.......situation 2 |
| call r_init(p0,nxyz,0.d0) |
| do i=1,lx1 |
| p0(i,1,1)=tcpre(i,1) |
| p0(1,i,1)=tcpre(i,1) |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_e(2,pcmor_cor(2,element_size),temp) |
| |
| c.......situation 3 |
| call r_init(p0,nxyz,0.d0) |
| do i=1,lx1 |
| p0(i,1,1)=tcpre(i,1) |
| p0(1,i,1)=tcpre(i,1) |
| p0(1,1,i)=tcpre(i,1) |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_e(3,pcmor_cor(3,element_size),temp) |
| |
| c.......situation 4 |
| call r_init(p0,nxyz,0.d0) |
| do j=1,lx1 |
| do i=1,lx1 |
| p0(i,j,1)=tcpre(i,j) |
| end do |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_f(4,pcmor_cor(4,element_size),temp) |
| |
| c.......situation 5 |
| call r_init(p0,nxyz,0.d0) |
| do j=1,lx1 |
| do i=1,lx1 |
| p0(i,j,1)=tcpre(i,j) |
| end do |
| end do |
| do i=1,lx1 |
| p0(1,1,i)=tcpre(i,1) |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_f(5,pcmor_cor(5,element_size),temp) |
| |
| c.......situation 6 |
| call r_init(p0,nxyz,0.d0) |
| do j=1,lx1 |
| do i=1,lx1 |
| p0(i,j,1)=tcpre(i,j) |
| p0(i,1,j)=tcpre(i,j) |
| end do |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_f(6,pcmor_cor(6,element_size),temp) |
| |
| c.......situation 7 |
| do j=1,lx1 |
| do i=1,lx1 |
| p0(i,j,1)=tcpre(i,j) |
| p0(i,1,j)=tcpre(i,j) |
| p0(1,i,j)=tcpre(i,j) |
| end do |
| end do |
| call laplacian(temp,p0,element_size) |
| call transfb_cor_f(7,pcmor_cor(7,element_size),temp) |
| |
| end do |
| c$OMP END PARALLEL DO |
| return |
| end |
| |
| |
| c------------------------------------------------------------------------ |
| subroutine setpcmo |
| c------------------------------------------------------------------------ |
| c compute the preconditioner by identifying its geometry configuration |
| c and sum the values from the precomputed elemental contributions |
| c------------------------------------------------------------------------ |
| |
| include 'header.h' |
| |
| integer face2, nb1, nb2, sizei, imor, enum, i,j, |
| & iel, iside, nn1, nn2 |
| |
| c$OMP PARALLEL DEFAULT(SHARED) PRIVATE(IMOR,IEL,ISIDE,I) |
| c$OMP DO |
| do imor=1,nvertex |
| ifpcmor(imor)=.false. |
| end do |
| c$OMP END DO nowait |
| |
| c$OMP DO |
| do iel=1,nelt |
| do iside=1,nsides |
| do i=1,4 |
| edgevis(i,iside,iel)=.false. |
| end do |
| end do |
| end do |
| c$OMP END DO |
| c$OMP END PARALLEL |
| |
| c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(IEL,iside,sizei, |
| c$OMP& imor,enum,face2,nb1,nb2,i,j,nn1,nn2) |
| |
| do iel=1,nelt |
| do iside=1,nsides |
| c.........for nonconforming faces |
| if(cbc(iside,iel).eq.3)then |
| sizei=size_e(iel) |
| |
| c...........vertices |
| |
| c...........ifpcmor(imor)=.true. indicates that mortar point imor has |
| c been visited |
| imor=idmo(1,1,1,1,iside,iel) |
| if(.not.ifpcmor(imor))then |
| c.............compute the preconditioner on mortar point imor |
| call pc_corner(imor) |
| ifpcmor(imor)=.true. |
| end if |
| |
| imor=idmo(lx1,1,1,2,iside,iel) |
| if(.not.ifpcmor(imor))then |
| call pc_corner(imor) |
| ifpcmor(imor)=.true. |
| end if |
| |
| imor=idmo(1,lx1,2,1,iside,iel) |
| if(.not.ifpcmor(imor))then |
| call pc_corner(imor) |
| ifpcmor(imor)=.true. |
| end if |
| |
| imor=idmo(lx1,lx1,2,2,iside,iel) |
| if(.not.ifpcmor(imor))then |
| call pc_corner(imor) |
| ifpcmor(imor)=.true. |
| end if |
| |
| c...........edges on nonconforming faces, enum is local edge number |
| do enum=1,4 |
| |
| c.............edgevis(enum,iside,iel)=.true. indicates that local edge |
| c enum of face iside of iel has been visited |
| if(.not.edgevis(enum,iside,iel))then |
| edgevis(enum,iside,iel)=.true. |
| |
| c...............Examing neighbor element information, |
| c calculateing the preconditioner value. |
| face2= f_e_ef(enum,iside) |
| if(cbc(face2,iel).eq.2)then |
| nb1=sje(1,1,face2,iel) |
| if(cbc(iside,nb1).eq.2)then |
| |
| c...................Compute the preconditioner on local edge enum on face |
| c iside of element iel, 1 is neighborhood information got |
| c by examing neighbors(nb1). For detailed meaning of 1, |
| c see subroutine com_dpc. |
| |
| call com_dpc(iside,iel,enum,1,sizei) |
| nb2=sje(1,1,iside,nb1) |
| edgevis(op(e_face2(enum,iside)), |
| & jjface(face2),nb2)=.true. |
| |
| elseif(cbc(iside,nb1).eq.3)then |
| call com_dpc(iside,iel,enum,2,sizei) |
| edgevis(op(enum),iside,nb1)=.true. |
| end if |
| |
| elseif(cbc(face2,iel).eq.3)then |
| edgevis(e_face2(enum,iside),face2,iel)=.true. |
| nb1=sje(1,2,face2,iel) |
| if(cbc(iside,nb1).eq.1)then |
| call com_dpc(iside,iel,enum,3,sizei) |
| nb2=sje(1,1,iside,nb1) |
| edgevis(op(enum),jjface(iside),nb2)=.true. |
| edgevis(op(e_face2(enum,iside)), |
| & jjface(face2),nb2)=.true. |
| elseif(cbc(iside,nb1).eq.2)then |
| call com_dpc(iside,iel,enum,4,sizei) |
| end if |
| else if (cbc(face2,iel).eq.0)then |
| call com_dpc(iside,iel,enum,0,sizei) |
| end if |
| end if |
| end do |
| |
| c...........mortar element interior (not edge of mortar) |
| |
| do nn1=1,2 |
| do nn2=1,2 |
| do j=2,lx1-1 |
| do i=2,lx1-1 |
| imor=idmo(i,j,nn1,nn2,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,nn1,nn2,sizei)+ |
| & pcmor_c(i,j,sizei+1)) |
| end do |
| end do |
| end do |
| end do |
| |
| c...........for i,j=lx1 there are duplicated mortar points, so |
| c pcmor_c needs to be doubled or quadrupled |
| i=lx1 |
| do j=2,lx1-1 |
| imor=idmo(i,j,1,1,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,1,1,sizei)+ |
| & pcmor_c(i,j,sizei+1)*2.d0) |
| imor=idmo(i,j,2,1,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,2,1,sizei)+ |
| & pcmor_c(i,j,sizei+1)*2.d0) |
| end do |
| |
| j=lx1 |
| imor=idmo(i,j,1,1,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,1,1,sizei)+ |
| & pcmor_c(i,j,sizei+1)*4.d0) |
| do i=2,lx1-1 |
| imor=idmo(i,j,1,1,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,1,1,sizei)+ |
| & pcmor_c(i,j,sizei+1)*2.d0) |
| imor=idmo(i,j,1,2,iside,iel) |
| dpcmor(imor) = 1.d0/(pcmor_nc1(i,j,1,2,sizei)+ |
| & pcmor_c(i,j,sizei+1)*2.d0) |
| end do |
| |
| end if |
| end do |
| end do |
| c$OMP END PARALLEL DO |
| |
| return |
| end |
| |
| c-------------------------------------------------------------------------- |
| subroutine pc_corner(imor) |
| c------------------------------------------------------------------------ |
| c calculate preconditioner value for vertex with mortar index imor |
| c------------------------------------------------------------------------ |
| |
| include 'header.h' |
| |
| double precision tmortemp |
| integer imor, inemo,ie, sizei,cornernumber, |
| & sface,sedge,iiface,iface,iiedge,iedge,n |
| |
| tmortemp=0.d0 |
| c.....loop over all elements sharing this vertex |
| do inemo=1,nemo(imor) |
| ie=emo(1,inemo,imor) |
| sizei=size_e(ie) |
| cornernumber=emo(2,inemo,imor) |
| sface=0 |
| sedge=0 |
| do iiface=1,3 |
| iface=f_c(iiface,cornernumber) |
| c.........sface sums the number of nonconforming faces sharing this vertex on |
| c one element |
| if(cbc(iface,ie).eq.3)then |
| sface=sface+1 |
| end if |
| end do |
| c.......sedge sums the number of nonconforming edges sharing this vertex on |
| c one element |
| do iiedge=1,3 |
| iedge=e_c(iiedge,cornernumber) |
| if(ncon_edge(iedge,ie))sedge=sedge+1 |
| end do |
| |
| c.......each n indicates how many nonconforming faces and nonconforming |
| c edges share this vertex on an element, |
| |
| if(sface.eq.0)then |
| if(sedge.eq.0)then |
| n=8 |
| elseif(sedge.eq.1)then |
| n=1 |
| elseif(sedge.eq.2)then |
| n=2 |
| elseif(sedge.eq.3)then |
| n=3 |
| end if |
| elseif (sface.eq.1)then |
| if (sedge.eq.1)then |
| n=5 |
| else |
| n=4 |
| end if |
| else if (sface.eq.2)then |
| n=6 |
| else if(sface.eq.3)then |
| n=7 |
| end if |
| |
| c.......sum the intermediate pre-computed preconditioner values for |
| c all elements |
| tmortemp=tmortemp+pcmor_cor(n,sizei) |
| |
| end do |
| |
| c.....dpcmor(imor) is the value of the preconditioner on mortar point imor |
| dpcmor(imor)=1.d0/tmortemp |
| |
| return |
| end |
| |
| c------------------------------------------------------------------------ |
| subroutine com_dpc(iside,iel,enumber,n,isize) |
| c------------------------------------------------------------------------ |
| c Compute preconditioner for local edge enumber of face iside |
| c on element iel. |
| c isize is element size, |
| c n is one of five different configurations |
| c anc1, ac, anc2, anc0 are coefficients for different edges. |
| c nc0 refers to nonconforming edge shared by two conforming faces |
| c nc1 refers to nonconforming edge shared by one nonconforming face |
| c nc2 refers to nonconforming edges shared by two nonconforming faces |
| c c refers to conforming edge |
| c------------------------------------------------------------------------ |
| |
| include 'header.h' |
| |
| integer n, isize,iside,iel, enumber, nn1start, nn1end, nn2start, |
| & nn2end, jstart, jend, istart, iend, i, j, nn1, nn2, imor |
| double precision anc1,ac,anc2,anc0,temp |
| |
| c.....different local edges have different loop ranges |
| if(enumber.eq.1)then |
| nn1start=1 |
| nn1end=1 |
| nn2start=1 |
| nn2end=2 |
| jstart=1 |
| jend=1 |
| istart=2 |
| iend=lx1-1 |
| elseif (enumber.eq.2) then |
| nn1start=1 |
| nn1end=2 |
| nn2start=2 |
| nn2end=2 |
| jstart=2 |
| jend=lx1-1 |
| istart=lx1 |
| iend=lx1 |
| elseif (enumber.eq.3) then |
| nn1start=2 |
| nn1end=2 |
| nn2start=1 |
| nn2end=2 |
| jstart=lx1 |
| jend=lx1 |
| istart=2 |
| iend=lx1-1 |
| elseif (enumber.eq.4) then |
| nn1start=1 |
| nn1end=2 |
| nn2start=1 |
| nn2end=1 |
| jstart=2 |
| jend=lx1-1 |
| istart=1 |
| iend=1 |
| end if |
| |
| c.....among the four elements sharing this edge |
| |
| c.....one has a smaller size |
| if(n.eq.1)then |
| anc1=2.d0 |
| ac=1.d0 |
| anc0=1.d0 |
| anc2=0.d0 |
| |
| c.....two (neighbored by a face) are of smaller size |
| else if (n.eq.2)then |
| anc1=2.d0 |
| ac=2.d0 |
| anc0=0.d0 |
| anc2=0.d0 |
| |
| c.....two (neighbored by an edge) are of smaller size |
| else if (n.eq.3)then |
| anc2=2.d0 |
| ac=2.d0 |
| anc1=0.d0 |
| anc0=0.d0 |
| |
| c.....three are of smaller size |
| else if (n.eq.4)then |
| anc1=0.d0 |
| ac=3.d0 |
| anc2=1.d0 |
| anc0=0.d0 |
| |
| c.....on the boundary |
| else if (n.eq.0)then |
| anc1=1.d0 |
| ac=1.d0 |
| anc2=0.d0 |
| anc0=0.d0 |
| end if |
| |
| c.....edge interior |
| do nn2=nn2start,nn2end |
| do nn1=nn1start,nn1end |
| do j=jstart,jend |
| do i=istart,iend |
| imor=idmo(i,j,nn1,nn2,iside,iel) |
| temp=anc1* pcmor_nc1(i,j,nn1,nn2,isize) + |
| & ac* pcmor_c(i,j,isize+1)+ |
| & anc0* pcmor_nc0(i,j,nn1,nn2,isize)+ |
| & anc2*pcmor_nc2(i,j,nn1,nn2,isize) |
| dpcmor(imor)=1.d0/temp |
| end do |
| end do |
| end do |
| end do |
| |
| c.......local edge 1 |
| if (enumber.eq.1) then |
| imor=idmo(lx1,1,1,1,iside,iel) |
| temp=anc1* pcmor_nc1(lx1,1,1,1,isize) + |
| & ac* pcmor_c(lx1,1,isize+1)*2.d0+ |
| & anc0* pcmor_nc0(lx1,1,1,1,isize)+ |
| & anc2*pcmor_nc2(lx1,1,1,1,isize) |
| c.......local edge 2 |
| elseif (enumber.eq.2) then |
| imor=idmo(lx1,lx1,1,2,iside,iel) |
| temp=anc1* pcmor_nc1(lx1,lx1,1,2,isize) + |
| & ac* pcmor_c(lx1,lx1,isize+1)*2.d0+ |
| & anc0* pcmor_nc0(lx1,lx1,1,2,isize)+ |
| & anc2*pcmor_nc2(lx1,lx1,1,2,isize) |
| c.......local edge 3 |
| elseif (enumber.eq.3) then |
| imor=idmo(lx1,lx1,2,1,iside,iel) |
| temp=anc1* pcmor_nc1(lx1,lx1,2,1,isize) + |
| & ac* pcmor_c(lx1,lx1,isize+1)*2.d0+ |
| & anc0* pcmor_nc0(lx1,lx1,2,1,isize)+ |
| & anc2*pcmor_nc2(lx1,lx1,2,1,isize) |
| c.......local edge 4 |
| elseif (enumber.eq.4) then |
| imor=idmo(1,lx1,1,1,iside,iel) |
| temp=anc1* pcmor_nc1(1,lx1,1,1,isize) + |
| & ac* pcmor_c(1,lx1,isize+1)*2.d0+ |
| & anc0* pcmor_nc0(1,lx1,1,1,isize)+ |
| & anc2*pcmor_nc2(1,lx1,1,1,isize) |
| end if |
| |
| dpcmor(imor)=1.d0/temp |
| |
| return |
| end |
| |
