@fem2d.pro

s3 = sqrt(3.)/2
x  = [0,1,2,3,4,0.5,1.5,2.5,3.5]
y  = [0,0,0,0,0,s3,s3,s3,s3]
x0 = x
y0 = y
;ixfix = [0]
ixfix = [0,4]
iyfix = [0,4]
ii1 = [0,1,2,3,0,5,1,6,2,7,3,8,5,6,7]
ii2 = [1,2,3,4,5,1,6,2,7,3,8,4,6,7,8]
nnode = n_elements(x)
nel   = n_elements(ii1)

r   = 0.01 + dblarr(nel)
e   = 100. + dblarr(nel)
rho = 1. + dblarr(nel)
g   = 1.
x0  = x
y0  = y


q   = fem2d(x0,y0,x,y,ixfix,iyfix,ii1,ii2,r,e,rho,g,$
             matrix=matrix,rhs=rhs,result=result)
x1  = x0+q.dx
y1  = y0+q.dy

window,0,xsize=600,ysize=600
position=[0.1,0.95,0.1,0.95]
plot,x,y,psym=6,yrange=[-2,3],xrange=[-0.5,4.5],/xs
for i=0,nel-1 do oplot,[x0[ii1[i]],x0[ii2[i]]],[y0[ii1[i]],y0[ii2[i]]]
oplot,x+q.dx,y+q.dy,psym=4
for i=0,nel-1 do oplot,[x1[ii1[i]],x1[ii2[i]]],[y1[ii1[i]],y1[ii2[i]]],line=2

;print,q.dy[2]
;for k=1,5 do begin
;    q2  = fem2d(x0,y0,x1,y1,ixfix,iyfix,ii1,ii2,r,e,rho,g,$
;             matrix=matrix,rhs=rhs,result=result)
;    print,q2.dy[2]
;    x1  = x0+q2.dx
;    y1  = y0+q2.dy
;    oplot,x1,y1,psym=4
;    for i=0,nel-1 do oplot,[x1[ii1[i]],x1[ii2[i]]],[y1[ii1[i]],y1[ii2[i]]],line=2
;endfor
end