module defparam_C implicit none integer,parameter::KTE = 20000 integer,parameter::KTJ = 20000 integer,parameter::KED = 20000 integer,parameter::KCM = 100 end module defparam_C !----------------------------------------------------------------------------------------- program f90_FEM_MESH4C !----------------------------------------------------------------------------------------- use defparam_C implicit none integer::NODE =0 integer::NELM =0 integer::mtj(1:KTE,1:4) =0 real(8)::px(1:KTJ) =0.0D0 real(8)::py(1:KTJ) =0.0D0 integer::ifix(1:KTJ) =0 character(len=50)::fnameR character(len=50)::fnameW call getarg(1,fnameR) call getarg(2,fnameW) call SQINPUT(NODE, NELM, mtj, px, py, ifix, fnameR) call SQLAPLAS(NODE, NELM, mtj, px, py, ifix) call SQDATA(NODE, NELM, mtj, px, py, fnameW) stop contains !----------------------------------------------------------------------------------------- subroutine SQINPUT(NODE, NELM, mtj, px, py, ifix, fnameR) !----------------------------------------------------------------------------------------- integer,intent(inout)::NODE integer,intent(inout)::NELM integer,intent(inout)::mtj(1:KTE,1:4) real(8),intent(inout)::px(1:KTJ) real(8),intent(inout)::py(1:KTJ) integer,intent(inout)::ifix(1:KTJ) character(len=50),intent(in)::fnameR integer::i,j open(11,file=fnameR,status='old') read(11,*) NODE,NELM do i = 1,NELM read(11,*) (mtj(i,j),j=1,4) end do do i = 1,NODE read(11,*) px(i),py(i),ifix(i) end do close(11) end subroutine SQINPUT !----------------------------------------------------------------------------------------- subroutine SQLAPLAS(NODE, NELM, mtj, px, py, ifix) !----------------------------------------------------------------------------------------- integer,intent(in)::NODE integer,intent(in)::NELM integer,intent(inout)::mtj(1:KTE,1:4) real(8),intent(inout)::px(1:KTJ) real(8),intent(inout)::py(1:KTJ) integer,intent(inout)::ifix(1:KTJ) integer::i,j,k,l,ja,jb,ip,ix,ka,kb,la,lb,itera,it,j1,j2,j3,j4,ielm integer::iflag,iarea real(8)::gx,gy,ar,s,xc,yc,cgrax,cgray integer::jac(1:KED,1:4) =0 integer::ihen(1:KED,1:2) =0 integer::jhen(1:KED) =0 integer::khen(1:KED) =0 integer::jnb(1:KED) =0 integer::nei(1:KED,1:KCM)=0 ix=0 do j = 1,NELM do k = 1,4 ka = mtj(j, k) kb = mtj(j, mod(k,4) + 1) iflag=0 ix = ix + 1 do l = 1,ix la = ihen(l, 1) lb = ihen(l, 2) if((ka == lb).and.(kb == la))then iflag=1 exit end if end do if(iflag==0)then ihen(ix, 1) = ka ihen(ix, 2) = kb jhen(ix) = j khen(ix) = k else jac(j, k) = jhen(l) jac(jhen(l), khen(l)) = j end if end do end do do i = 1,NELM do j = 1,4 if(jac(i, j) == 0)then ja = mtj(i, j) jb = mtj(i, mod(j,4) + 1) ifix(ja) = 1 ifix(jb) = 1 end if end do end do do i = 1,NELM do j = 1,4 ip = mtj(i, j) jnb(ip) = jnb(ip) + 1 nei(ip, jnb(ip)) = i end do end do itera = 5 do it = 1,itera do i = 1,NODE if(ifix(i) == 0)then gx = 0.0D0 gy = 0.0D0 ar = 0.0D0 do j = 1,jnb(i) ielm = nei(i, j) j1 = mtj(ielm, 1) j2 = mtj(ielm, 2) j3 = mtj(ielm, 3) j4 = mtj(ielm, 4) s = SQAREA(ielm, mtj, px, py, iarea) xc = (px(j1) + px(j2) + px(j3) + px(j4)) / 4.0 yc = (py(j1) + py(j2) + py(j3) + py(j4)) / 4.0 ar = ar + s gx = gx + s * xc gy = gy + s * yc end do cgrax = gx / ar cgray = gy / ar px(i) = cgrax py(i) = cgray end if end do end do end subroutine SQLAPLAS !----------------------------------------------------------------------------------------- subroutine SQDATA(NODE, NELM, mtj, px, py, fnameW) !----------------------------------------------------------------------------------------- integer,intent(in)::NODE integer,intent(in)::NELM integer,intent(in)::mtj(1:KTE,1:4) real(8),intent(in)::px(1:KTJ) real(8),intent(in)::py(1:KTJ) character(len=50),intent(in)::fnameW integer::i,j,iarea real(8)::area open(12, file=fnameW, status='replace') write(12,'(2(i5))') NODE,NELM do i = 1,NELM write(12,'(7(i5))') (mtj(i,j),j=1,4) end do do i = 1,NODE write(12,'(2(e15.7))') px(i),py(i) end do ! Find irregular elements do i=1,NELM area=SQAREA(i, mtj, px, py, iarea) if(iarea==1)then write(12,'("Irregular element: No,area,iarea=",i5,", ",e15.7,", ",i5)'),i,area,iarea do j=1,4 write(12,'(i5,2(e15.7))'),mtj(i,j),px(mtj(i,j)),py(mtj(i,j)) end do end if end do close(12) end subroutine SQDATA !----------------------------------------------------------------------------------------- real(8) function SQAREA(ielm, mtj, px, py, iarea) !----------------------------------------------------------------------------------------- integer,intent(in)::ielm integer,intent(in)::mtj(1:KTE,1:4) real(8),intent(in)::px(1:KTJ) real(8),intent(in)::py(1:KTJ) integer,intent(out)::iarea integer::j1,j2,j3,j4 real(8)::area1,area2,area3,area4 iarea=1 j1 = mtj(ielm, 1) j2 = mtj(ielm, 2) j3 = mtj(ielm, 3) j4 = mtj(ielm, 4) area1 = 0.5D0 * & & (px(j1) * py(j2) + px(j2) * py(j3) + px(j3) * py(j1) - px(j1) * py(j3) - px(j2) * py(j1) - px(j3) * py(j2)) area2 = 0.5D0 * & & (px(j1) * py(j3) + px(j3) * py(j4) + px(j4) * py(j1) - px(j1) * py(j4) - px(j3) * py(j1) - px(j4) * py(j3)) j1 = mtj(ielm, 2) j2 = mtj(ielm, 3) j3 = mtj(ielm, 4) j4 = mtj(ielm, 1) area3 = 0.5D0 * & & (px(j1) * py(j2) + px(j2) * py(j3) + px(j3) * py(j1) - px(j1) * py(j3) - px(j2) * py(j1) - px(j3) * py(j2)) area4 = 0.5D0 * & & (px(j1) * py(j3) + px(j3) * py(j4) + px(j4) * py(j1) - px(j1) * py(j4) - px(j3) * py(j1) - px(j4) * py(j3)) if((0.0D0