currentdir(); mkdir( "repos" ); currentdir("repos"); currentdir(); march( 'create', ".", 200 ); Plotter:=module() description "Graphic Package for Plotter Simulation"; export start_2D_plot,start_3D_plot,start_2D_animation,start_3D_animation,\ show_2D_plot,show_3D_plot,show_2D_animation,show_3D_animation,\ upgrade_frames,n_arrow,new_arrow,pen_up,pen_down,move_pen,plot_element,\ new_axes,new_color,new_thickness,new_pen_style,new_plot_style,\ new_plot_shading,draw_symbol,draw_text,draw_title,draw_label,new_scaling; local view_plot; option package; start_2D_plot:=proc() global x_view_max,y_view_max,x_view_min,y_view_min,plot_dimension,xyz_coor,\ pen_color,pen_thickness,pen_position,pen_style,plot_style,plot_title,\ plot_axes,plot_scaling,draw; x_view_max,y_view_max:=-1000,-1000; x_view_min,y_view_min:=1000,1000; plot_dimension:=2; xyz_coor:=[0,0]; pen_color:=COLOR(RGB,0,0,0); pen_thickness:=THICKNESS(0); pen_style:=LINESTYLE(0); pen_position:=0; plot_style:=STYLE(PATCH); plot_title:=NULL; plot_axes:=AXESSTYLE(NONE); plot_scaling:=SCALING(CONSTRAINED); draw:={}; NULL end proc: start_3D_plot:=proc() global x_view_max,y_view_max,x_view_min,y_view_min,plot_dimension,xyz_coor,\ pen_position,pen_style,pen_color,pen_thickness,plot_style,plot_title,\ plot_shading,plot_axes,draw; x_view_max,y_view_max:=-1000,-1000; x_view_min,y_view_min:=1000,1000; plot_dimension:=3; xyz_coor:=[0,0,0]; pen_position:=0; pen_style:=LINESTYLE(0); pen_color:=COLOR(ZSHADING); pen_thickness:=THICKNESS(0); plot_style:=STYLE(PATCH); plot_shading:=SHADING(NONE); plot_title:=NULL; plot_axes:=AXESSTYLE(NONE); draw:={}; NULL end proc: start_2D_animation:=proc() global x_view_max,y_view_max,x_view_min,y_view_min,plot_dimension,xyz_coor,\ pen_position,pen_style,pen_color,pen_thickness,plot_style,plot_title,\ plot_axes,plot_scaling,draw,frames; x_view_max,y_view_max:=-1000,-1000; x_view_min,y_view_min:=1000,1000; plot_dimension:=2; xyz_coor:=[0,0]; pen_position:=0; pen_style:=LINESTYLE(0); pen_color:=COLOR(RGB,0,0,0); pen_thickness:=THICKNESS(0); plot_style:=STYLE(PATCH); plot_title:=NULL; plot_axes:=AXESSTYLE(NONE); plot_scaling:=SCALING(CONSTRAINED); draw:={}; frames:=[]; NULL end proc: start_3D_animation:=proc() global x_view_max,y_view_max,x_view_min,y_view_min,plot_dimension,xyz_coor,\ pen_position,pen_style,pen_color,pen_thickness,plot_style,plot_title,\ plot_shading,plot_axes,draw,frames; x_view_max,y_view_max:=-1000,-1000; x_view_min,y_view_min:=1000,1000; plot_dimension:=3; xyz_coor:=[0,0,0]; pen_position:=0; pen_style:=LINESTYLE(0); pen_color:=COLOR(ZSHADING); pen_thickness:=THICKNESS(0); plot_style:=STYLE(PATCH); plot_shading:=SHADING(NONE); plot_title:=NULL; plot_axes:=AXESSTYLE(NONE); draw:={}; frames:=[]; NULL end proc: show_2D_plot:=proc(llc::list(algebraic),urc::list(algebraic)) global draw,plot_title,plot_axes,plot_scaling; local opts; opts:=plot_title,plot_axes,plot_scaling; if nargs=2 then opts:=opts,VIEW(llc[1]..urc[1],llc[2]..urc[2]) else opts:=opts,VIEW(DEFAULT,DEFAULT) end if; PLOT(op(draw),opts) end proc: show_3D_plot:=proc(llc::list(algebraic),urc::list(algebraic)) global draw,plot_shading,pen_color,plot_axes,plot_title; local opts; opts:=plot_title,plot_axes; if plot_shading<>SHADING(NONE) then opts:=opts,ORIENTATION(25.,60.); if pen_color<>COLOR(ZSHADING) then opts:=opts,LIGHTMODEL(LIGHT_1) end if end if; if nargs=2 then opts:=opts,VIEW(llc[1]..urc[1],llc[2]..urc[2],DEFAULT) else opts:=opts,VIEW(DEFAULT,DEFAULT,DEFAULT) end if; PLOT3D(op(draw),opts) end proc: show_2D_animation:=proc(llc::list(algebraic),urc::list(algebraic)) global frames,plot_title,plot_axes,plot_scaling; local opts; opts:=plot_axes,plot_title,plot_scaling; if nargs=2 then opts:=opts,VIEW(llc[1]..urc[1],llc[2]..urc[2]) else opts:=opts,VIEW(DEFAULT,DEFAULT) end if; PLOT(ANIMATE(frames[]),opts) end proc: show_3D_animation:=proc(llc::list(algebraic),urc::list(algebraic)) global frames,pen_color,plot_shading,plot_title,plot_axes; local opts; opts:=plot_axes,plot_title; if plot_shading<>SHADING(NONE) then opts:=opts,ORIENTATION(25.,60.); if pen_color<>COLOR(ZSHADING) then opts:=opts,LIGHTMODEL(LIGHT_1) end if end if; if nargs=2 then opts:=opts,VIEW(llc[1]..urc[1],llc[2]..urc[2],DEFAULT) else opts:=opts,VIEW(DEFAULT,DEFAULT,DEFAULT) end if; PLOT3D(ANIMATE(frames[]),opts) end proc: upgrade_frames:=proc() global frames,draw; frames:=[frames[],[op(draw)]]; draw:={}; NULL end proc: n_arrow:=proc(arg::anything) global pen_color,pen_thickness,pen_style,draw; local opts,i,j,point,vector,aopts; opts:=pen_color,pen_thickness,pen_style; aopts:={}; for i in args[] do j:=convert(op(1,i),string); if j="point" then point:=op(2,i) elif j="vector" then vector:=op(2,i) else aopts:={op(aopts),i} end if end do; view_plot(point[1]+vector[1],point[2]+vector[2]); op(plots[arrow](point,vector,op(aopts))): draw:={op(draw),%}; NULL end proc: new_arrow:=proc(point::list,vect::list,t1::algebraic,t2::algebraic,t3::algebraic) # [point], [vect], a_base_width, a_head_width, r_head_length global pen_color,pen_thickness,pen_style,draw; local opts,x,y,a,b,L,Cos,Sin,v,i; opts:=pen_color,pen_thickness,pen_style; x:=point[1];y:=point[2];a:=vect[1];b:=vect[2]; view_plot(x+a,y+b); if has(vect,'undefined') or (abs(a)<10^(1-Digits) and abs(b)<10^(1-Digits)) then POLYGONS([]) else L:=evalf(sqrt(a^2+b^2)); Cos:=evalf(a/L);Sin:=evalf(b/L); v[1]:=[x+t1*Sin/2,y-t1*Cos/2];v[2]:=[x-t1*Sin/2,y+t1*Cos/2]; v[3]:=[x-t1*Sin/2-t3*Cos*L+a,y+t1*Cos/2-t3*Sin*L+b]; v[4]:=[x+t1*Sin/2-t3*Cos*L+a,y-t1*Cos/2-t3*Sin*L+b]; v[5]:=[x-t2*Sin/2-t3*Cos*L+a,y+t2*Cos/2-t3*Sin*L+b]; v[6]:=[x+a,y+b];v[7]:=[x+t2*Sin/2-t3*Cos*L+a,y-t2*Cos/2-t3*Sin*L+b]; v:=seq(evalf(v[i]),i=1..7); draw:={op(draw),POLYGONS([v[1],v[2],v[3],v[4]],[v[5],v[6],v[7]],opts[]),\ CURVES([v[1],v[2],v[3],v[5],v[6],v[7],v[4],v[1]])}; NULL end if end proc: view_plot:=proc(x::algebraic,y::algebraic) global x_view_max,y_view_max,x_view_min,y_view_min; x_view_max:=max(x_view_max,x); y_view_max:=max(y_view_max,y); x_view_min:=min(x_view_min,x); y_view_min:=min(y_view_min,y); NULL end proc: pen_up:=proc() global pen_position; pen_position:=0; NULL end proc: pen_down:=proc() global pen_position; pen_position:=1; NULL end proc: move_pen:=proc(p::list(algebraic),pen_pos::nonnegint) global xyz_coor,pen_color,pen_thickness,pen_position,pen_style,draw; local opts,xyz,i; opts:=pen_color,pen_thickness,pen_style; xyz:=xyz_coor; xyz_coor:=map(i->i,p); if nargs=2 then pen_position:=min(pen_pos,1) fi; view_plot(xyz_coor[1],xyz_coor[2]); if pen_position=1 then draw:={CURVES([[map(i->evalf(i),xyz)[]],[map(i->evalf(i),xyz_coor)[]]],\ opts),op(draw)} end if; NULL end proc: plot_element:=proc(x::list(algebraic),y::list(algebraic),z::list(algebraic)) global pen_color,pen_thickness,pen_style,draw; local vertices,opts,i; opts:=pen_style,pen_color,pen_thickness,plot_style; if nargs>2 then vertices:=seq([x[i],y[i],z[i]],i=1..nops(x)) else vertices:=seq([x[i],y[i]],i=1..nops(x)) end if; draw:={op(draw),POLYGONS([vertices],opts)}; NULL end proc: new_axes:=proc(ax::nonnegint) global plot_axes; local st; if nargs=1 then if ax=0 then st:=NONE elif ax=1 then st:=NORMAL else st:=BOX fi; else st:=NONE end if; plot_axes:=AXESSTYLE(st); NULL end proc: new_scaling:=proc(sc::anything) global plot_scaling; local st; if nargs=1 then st:=UNCONSTRAINED else st:=CONSTRAINED end if; plot_scaling:=SCALING(st); NULL end proc: new_color:=proc(r::algebraic,g::algebraic,b::algebraic) global plot_dimension,pen_color; local color; if nargs=3 then color:=[RGB,min(1,max(0,r)),min(1,max(0,g)),min(1,max(0,b))] else if plot_dimension=2 then color:=[RGB,0,0,0] else color:=[ZSHADING] end if end if; pen_color:=COLOR(color[]); NULL end proc: new_thickness:=proc(t::nonnegint) global pen_thickness; local thick; if nargs=1 then thick:=min(3,t) else thick:=0 end if; pen_thickness:=THICKNESS(thick); NULL end proc: new_pen_style:=proc(st::nonnegint) global pen_style; local style; if nargs=1 then style:=min(7,st) else style:=0 end if; pen_style:=LINESTYLE(style); NULL end proc: new_plot_style:=proc(st::nonnegint) global plot_style; local style; if nargs=1 then style:=min(3,st); if style=1 then style:=[PATCH] elif style=2 then style:=[CONTOUR] else style:=[HIDDEN] end if else style:=[PATCH] end if; plot_style:=STYLE(style[]); NULL end proc: new_plot_shading:=proc(st::nonnegint) global plot_shading; local style; if nargs=1 then style:=min(3,st); if style=1 then style:=[ZSHADING] elif style=2 then style:=[XYSHADING] else style:=[XYZSHADING] end if else style:=[NONE] end if; plot_shading:=SHADING(style[]); NULL end proc: draw_symbol:=proc(p::list(algebraic),s::nonnegint,z::nonnegint) global pen_color,draw; local ss,symb; if nargs<2 then symb:='DEFAULT' else ss:=min(5,s); if ss=0 then symb:='DEFAULT' elif ss=1 then symb:='BOX' elif ss=2 then symb:='CROSS' elif ss=3 then symb:='CIRCLE' elif ss=4 then symb:='POINT' else symb:='DIAMOND' end if end if; view_plot(p[1],p[2]); if nargs=3 then draw:={op(draw),POINTS([map(i->evalf(i),p)[]],SYMBOL(symb,z),pen_color)} else draw:={op(draw),POINTS([map(i->evalf(i),p)[]],SYMBOL(symb),pen_color)}; end if; NULL end proc: draw_text:=proc(p::list(algebraic),text::anything,\ align::list(algebraic),fsize::nonnegint) global pen_color,draw; local size,h_align,v_align,align_opts; if nargs>3 then size:=fsize else size:=10 end if; view_plot(p[1],p[2]);h_align,v_align,align_opts:=0,0,[]; if nargs>2 then if align[1]<>0 then h_align:=min(1,abs(align[1]))*signum(align[1]) end if; if nops(align)=2 and align[2]<>0 then v_align:=min(1,abs(align[2]))*signum(align[2]) end if; if h_align=1 then h_align:='ALIGNRIGHT' elif h_align=-1 then h_align:='ALIGNLEFT' end if; if v_align=1 then v_align:='ALIGNABOVE' elif v_align=-1 then v_align:='ALIGNBELOW' end if; if h_align<>0 then align_opts:=[align_opts[],h_align] end if; if v_align<>0 then align_opts:=[align_opts[],v_align] end if end if; if align_opts<>[] then draw:={op(draw),TEXT([map(i->evalf(i),p)[]],\ 'text',align_opts[],pen_color,FONT(TIMES,BOLD,size))} else draw:={op(draw),TEXT([map(i->evalf(i),p)[]],\ 'text',pen_color,FONT(TIMES,BOLD,size))} end if; NULL end proc: draw_title:=proc(text::anything,fsize::nonnegint) global draw,plot_title; local size; if nargs>1 then size:=fsize else size:=14 end if; plot_title:=TITLE('text',FONT(TIMES,BOLD,size)); NULL end proc: draw_label:=proc(label::list(anything),col::listlist) global x_view_max,y_view_max,x_view_min,y_view_min; local dx,yy,lab; dx:=(x_view_max-x_view_min)/nops(label); yy:=y_view_min-(y_view_max-y_view_min)*0.1; for lab in label do member(lab,label,'i'); new_color(col[i][]); draw_text([x_view_min+(2*i-1)*dx/2,yy],lab,[0,0],12) end do; NULL end proc: end module: Cgt_fem:=module() description "2D Potential Flow Analysis with CGT Finite Elements"; export cgt_fem,read_save_data,compute_fluxes,element_geometry; local read_data,read_line,line_contents,save_data,\ initialize_global_matrices,assemble_global_matrices,\ assemble_cofficients,element_k,set_up_d_b,shape_f,\ lder_shape_f,jacob,gder_shape_f,assemble_point_sources,\ assemble_element_sources,assemble_boundary_velocities,\ force_exact_boundary_conditions,compute_potentials,\ check_boundary_fluxes,compute_velocities,element_velocities,\ element_fluxes,loop_123; option package; read_save_data:=proc() local reply; reply:=0; while reply<>y and reply<>n do reply:=readstat("read data from a file (y/n) ? ") end do; if reply=y then read_data() end if; reply:=0; while reply<>y and reply<>n do reply:=readstat("save data into a file (y/n) ? ") end do; if reply=y then save_data() end if; NULL end proc: read_data:=proc() global tcase,control,nods,elems,mat_props,bdr_conds,b_velts,p_srcs,e_srcs; local reply,fle,nod,elem,line,lne,i; reply:=readstat("file name: "); fle:=fopen(reply,READ,TEXT); control,nods,elems,mat_props,bdr_conds,b_velts,p_srcs,e_srcs:=seq([],i=1..8); lne:=readline(fle):line:=sscanf(lne,"%s"); while line<>["*end*"] do if line=["*control*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%s"); if substring(lne,1..5)="title" then tcase:=lne[6..length(lne)] elif substring(lne,1..13)="point sources" then control:=[control[],"point sources"] elif substring(lne,1..15)="element sources" then control:=[control[],"element sources"] elif substring(lne,1..19)="boundary velocities" then control:=[control[],"boundary velocities"] fi end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*nodes*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a"); if line<>0 and line<>[] then line_contents(line,[nods[],[seq(line[i],i=2..3)]],nods) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*elements*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a%a%a"); if line<>0 and line<>[] then line_contents(line,[elems[],[seq(line[i],i=2..5)]],elems) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*constraints*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a"); if line<>0 and line<>[] then line_contents(line,[bdr_conds[],[seq(line[i],i=1..2)]],bdr_conds) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif substring(lne,1..15)="*point sources*" then print(lne[1..15]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a"); if line<>0 and line<>[] then line_contents(line,[p_srcs[],[seq(line[i],i=1..2)]],p_srcs) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif substring(lne,1..21)="*boundary velocities*" then print(lne[1..21]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a"); if line<>0 and line<>[] then line_contents(line,[b_velts[],[seq(line[i],i=1..3)]],b_velts) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif substring(lne,1..17)="*element sources*" then print(lne[1..17]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a"); if line<>0 and line<>[] then line_contents(line,[e_srcs[],[seq(line[i],i=1..2)]],e_srcs) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*materials*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a%a"); if line<>0 and line<>[] then line_contents(line,[mat_props[],[seq(line[i],i=2..4)]],mat_props) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if else lne:=readline(fle):line:=sscanf(lne,"%s") end if end do; fclose(fle); NULL end proc: read_line:=proc(fle::integer,format::string) local lne; lne:=readline(fle); lne,sscanf(lne,format) end proc: line_contents:=proc(line::{list,integer},lst::list,var::evaln) if line<>0 and line<>[] then var:=lst else NULL end if end proc: save_data:=proc() global tcase,control,nods,elems,mat_props,bdr_conds,b_velts,p_srcs,e_srcs; local reply,file,i; reply:=readstat("file name: "); file:=fopen(reply,WRITE,TEXT); fprintf(file,"\n*control* [title,point sources,"); fprintf(file,"element sources,boundary velocities]\n"); fprintf(file,"%s%s\n","title",tcase); map(proc(i,file) fprintf(file,"%s\n",i) end,control,file); fprintf(file,"\n*materials* [material,x-permeability,y-permeability,"); fprintf(file,"angle of local x-direction]\n"); seq(fprintf(file,"%a %a %a %a\n",i,mat_props[i][]),i=1..nops(mat_props)); fprintf(file,"\n*nodes* [node,x,y]\n"); seq(fprintf(file,"%a %a %a\n",i,nods[i][]),i=1..nops(nods)); fprintf(file,"\n*elements* [element,node1,node2,node3,material]\n"); seq(fprintf(file,"%a %a %a %a %a\n",i,elems[i][]),i=1..nops(elems)); fprintf(file,"\n*constraints* [node,potential]\n"); map(proc(i,file) fprintf(file,"%a %a\n",i[]) end,bdr_conds,file); if nops(b_velts)>0 then fprintf(file,"\n*boundary velocities* [node1,node2,boundary-normal velocity]\n"); map(proc(i,file) fprintf(file,"%a %a %a\n",i[]) end,b_velts,file) end if; if nops(p_srcs)>0 then fprintf(file,"\n*point sources* [node,Q]\n"); map(proc(i,file) fprintf(file,"%a %a\n",i[]) end,p_srcs,file) end if; if nops(e_srcs)>0 then fprintf(file,"\n*element sources* [element,q]\n"); map(proc(i,file) fprintf(file,"%a %a\n",i[]) end,e_srcs,file) end if; fprintf(file,"\n*end*\n"); fclose(file); NULL end proc: cgt_fem:=proc() # *** main program *** global nods,elems,g_k,g_p,n_potls; local i,cpu_time,talk; cpu_time:=time(); if nargs=0 then talk:=NULL end if; initialize_global_matrices(talk); assemble_global_matrices(talk); force_exact_boundary_conditions(talk); n_potls:=linalg[linsolve](g_k,g_p); compute_potentials(talk); check_boundary_fluxes(talk); compute_velocities(talk); compute_fluxes([seq(i,i=1..nops(elems))],talk); cpu_time:=convert(time()-cpu_time,string); print(`Cpu time: `||cpu_time||` seconds`); NULL end proc: initialize_global_matrices:=proc() global nods,g_k,g_p,bvlts; local dim,i,j,p; dim:=nops(nods);bvlts:=[]; g_k:=array(symmetric,1..dim,1..dim,[seq([seq(0,j=1..i)],i=1..dim)]); g_p:=array([seq(0,j=1..dim)]); if nargs<>0 then print(`Initializing global matrices:`);print(g_k,g_p) end if; NULL end proc: assemble_global_matrices:=proc() global control; local c; assemble_cofficients(args); if control<>[] then for c in control do if c="point sources" then assemble_point_sources(args) elif c="element sources" then assemble_element_sources(args) elif c="boundary velocities" then assemble_boundary_velocities(args) end if end do end if; NULL end proc: assemble_cofficients:=proc() global nods,elems,mat_props,g_k; local e_k,ek,parms,e,n,mat,i,j,x,y,X,Y,Kx,Ky; e_k:=element_k([seq([X[i],Y[i]],i=1..3)],[Kx,Ky]); parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],Kx=mat[1],Ky=mat[2]]; for e in elems do mat:=mat_props[e[4]][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); ek:=simplify(subs(eval(parms)[],eval(e_k)),assume=positive); n:=e[1..3]; for i to 3 do for j to i do g_k[n[i],n[j]]:=g_k[n[i],n[j]]+ek[i,j] end do end do; if nargs<>0 then print('`Assembling element`'*e*'`:`'); print(simplify(g_k,assume=positive)) end if end do; NULL end proc: element_k:=proc(e_nods::list(list),e_props::list) local s,t,d,b,det_jac,kernel; d,b,det_jac:=set_up_d_b(e_nods,e_props,[s,t]); kernel:=linalg[multiply](linalg[transpose](b),d,b)*det_jac; evalm(int(int(kernel,t=0..1-s),s=0..1)) end proc: set_up_d_b:=proc(e_nods::list(list),e_props::list,nat_crds::list) local Kx,Ky,d,gder,s,t,b,det_jac,i; option remember; Kx,Ky:=e_props[1],e_props[2]; d:=array([[-Kx,0],[0,-Ky]]); gder:=gder_shape_f(e_nods,det_jac,[s,t]); b:=seq(array([[gder[i][1]],[gder[i][2]]]),i=1..3); b:=linalg[augment](b[1],b[2],b[3]); if nargs>2 then det_jac:=subs(s=nat_crds[1],t=nat_crds[2],det_jac); b:=subs(s=nat_crds[1],t=nat_crds[2],eval(b)) end if; map(eval,[d,b,det_jac])[] end proc: shape_f:=proc(nat_crds::list) local s,t,N; option remember; N:=[1-s-t,s,t]; # nodes at: 1->(0,0) 2->(1,0) 3->(0,1) if nargs<>0 then N:=subs(s=nat_crds[1],t=nat_crds[2],N) end if; eval(N) end proc: lder_shape_f:=proc(nat_crds::list) local N,s,t,lder,i; option remember; N:=shape_f([s,t]); lder:=[seq([diff(N[i],s),diff(N[i],t)],i=1..3)]; if nargs<>0 then lder:=subs(s=nat_crds[1],t=nat_crds[2],lder) end if; eval(lder) end proc: jacob:=proc(e_nods::list(list),nat_crds::list) local nod_crds,N,s,t,x,y,jac,i; option remember; nod_crds:=map(i->i[],e_nods); N:=shape_f([s,t]); N:=seq(array([[N[i],0],[0,N[i]]]),i=1..3); N:=linalg[augment](N[1],N[2],N[3]);linalg[multiply](N,nod_crds); x,y:=%[1],%[2]; jac:=linalg[transpose](linalg[jacobian]([x,y],[s,t])); if nargs>1 then jac:=subs(s=nat_crds[1],t=nat_crds[2],eval(jac)) end if; eval(jac); end proc: gder_shape_f:=proc(e_nods::list(list),det_jac::evaln,nat_crds::list) local jac,i_jac,s,t,gder,i; option remember; jac:=jacob(e_nods,[s,t]); det_jac:=linalg[det](jac); i_jac:=linalg[inverse](jac); gder:=[seq(linalg[multiply](i_jac,lder_shape_f([s,t])[i]),i=1..3)]; if nargs>2 then gder:=subs(s=nat_crds[1],t=nat_crds[2],eval(gder)) end if; eval(gder) end proc: assemble_point_sources:=proc() global p_srcs,g_p; local Q,i; for Q in p_srcs do i:=Q[1]; g_p[i]:=g_p[i]+Q[2] end do; if nargs<>0 then print(`Assembling point sources:`);print(g_p) end if; NULL end proc: assemble_element_sources:=proc() global nods,elems,e_srcs,g_p; local q,e,n,area,i,j,Q; for q in e_srcs do e:=q[1]; n:=[seq([seq(nods[elems[e][i]][j],j=1..2)],i=1..3)]; area:=element_geometry(n);Q:=q[2]*area/3; for i to 3 do j:=elems[e][i]; g_p[j]:=g_p[j]+Q end do end do; if nargs<>0 then print(`Assembling element sources:`);print(g_p) end if; NULL end proc: assemble_boundary_velocities:=proc() global b_velts,bvlts,nods,g_p; local b,m,n,q,i; bvlts:={}; for b in b_velts do m,n,q:=b[]; sqrt((nods[n][1]-nods[m][1])^2+(nods[n][2]-nods[m][2])^2); q:=q*%/2; g_p[m]:=g_p[m]+q; g_p[n]:=g_p[n]+q; bvlts:=bvlts union {m} union {n} end do; bvlts:=convert(bvlts,list); for i to nops(bvlts) do n:=bvlts[i];bvlts[i]:=[n,[g_p[n]]]; end do; if nargs<>0 then print(`Assembling boundary normal velocities:`);print(g_p) end if; NULL end proc: force_exact_boundary_conditions:=proc() global bflux,bdr_conds,nods,g_k,g_p; local dim,cond,n,i,j; dim:=nops(nods);bflux:=[]; for cond in bdr_conds do n:=cond[1]; bflux:=[bflux[],[n,[seq(g_k[n,j],j=1..dim)]]] end do; for cond in bdr_conds do n:=cond[1]; for i to dim do g_p[i]:=g_p[i]-cond[2]*g_k[i,n] end do; for i to dim do g_k[n,i]:=0 end do; g_k[n,n]:=1; g_p[n]:=cond[2] end do; if nargs<>0 then print(`Forcing exact boundary conditions:`); print(simplify(g_k,assume=positive)); print(simplify(g_p,assume=positive)) end if; NULL end proc: compute_potentials:=proc() global n_potls; n_potls:=simplify(convert(n_potls,list),assume=positive); if nargs<>0 then print(`Nodal potentials (n_potls):`);print(n_potls) end if; NULL end proc: check_boundary_fluxes:=proc() global bflux,bvlts,n_potls; local i,j,flx,tflx,r,p,Q,q,n,area; i,tflx:=0,0; for flx in bflux do i:=i+1; r:=add(flx[2][j]*n_potls[j],j=1..nops(n_potls)); bflux[i]:=[flx[1],[r]]; tflx:=tflx+r end do; for Q in p_srcs do tflx:=tflx+Q[2]; end do; for q in e_srcs do n:=[seq([seq(nods[elems[q[1]][i]][j],j=1..2)],i=1..3)]; area:=element_geometry(n);Q:=q[2]*area; tflx:=tflx+Q; end do; if bvlts<>[] then bflux:=[bflux[],bvlts[]]; for q in bvlts do tflx:=tflx+q[2][]; end do; end if; bflux:=simplify(bflux,assume=positive); tflx:=simplify(tflx,assume=positive); if nargs<>0 then print(`Boundary fluxes (bflux):`*['node',['flux']]); print(bflux[]);tflx:=convert(tflx,string); print(`Flux equilibrium: `*'Sigma*fluxes'*`= `||tflx) end if; NULL end proc: compute_velocities:=proc() global mat_props,nods,n_potls,e_grads,e_velts,n_velts,n_kined,total_internal_energy; local dim,sig,grd,kin,kine,sigma,parms,hit,e,mat,x,y,v,i,j,X,Y,V,Kx,Ky,Ang; element_velocities([seq([X[i],Y[i]],i=1..3)],[Kx,Ky,Ang],[seq(V[i],i=1..3)]); sig,grd,kin:=%; parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],Kx=mat[1],Ky=mat[2],Ang=mat[3]]; parms:=[parms[],seq(V[i]=v[i],i=1..3)]; dim:=nops(nods);hit:=[seq(0,i=1..dim)]; n_kined:=[seq(0,i=1..dim)];total_internal_energy:=0; e_grads:=[];e_velts:=[];n_velts:=[seq([0,0],i=1..dim)]; for e in elems do mat:=mat_props[e[4]][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); v:=[seq(n_potls[e[i]],i=1..3)]; subs(eval(parms)[],eval(grd)); e_grads:=[e_grads[],%]; sigma:=subs(eval(parms)[],eval(sig)); e_velts:=[e_velts[],sigma]; kine:=subs(eval(parms)[],eval(kin)); total_internal_energy:=total_internal_energy+kine[1]; for i in e[1..3] do n_velts[i]:=[seq(n_velts[i][j]+sigma[j],j=1..2)]; n_kined[i]:=n_kined[i]+kine[2]; hit[i]:=hit[i]+1 end do end do; n_velts:=[seq([seq(n_velts[i][j]/hit[i],j=1..2)],i=1..dim)]; n_kined:=[seq(n_kined[i]/hit[i],i=1..dim)]; n_kined:=simplify(n_kined,assume=positive); e_velts:=simplify(e_velts,assume=positive); n_velts:=simplify(n_velts,assume=positive); e_grads:=simplify(e_grads,assume=positive); total_internal_energy:=simplify(total_internal_energy,assume=positive); if nargs<>0 then sig:=['v[x]','v[y]']:grd:=['grad[x]','grad[y]']:kin:='w'; print(`Element gradients (e_grads): `*grd);print(e_grads[]); print(`Element velocities (e_velts): `*sig);print(e_velts[]); print(`Nodal velocities (n_velts): `*sig);print(n_velts[]); print(`Nodal internal energy density (n_kined): `*kin);print(n_kined[]); print(`Total internal energy: `);print(total_internal_energy) end if; NULL end proc: element_velocities:=proc(e_nods::list(list),e_props::list,u::list) local d,b,det_jac,lgrad,grad,lsigma,T,sn,cn,ang,sigma,w,U; d,b,det_jac:=set_up_d_b(e_nods,e_props); grad:=linalg[multiply](b,u); #global ang:=e_props[3]*Pi/180;sn:=sin(ang);cn:=cos(ang); T:=array([[cn,sn],[-sn,cn]]); lgrad:=linalg[multiply](linalg[transpose](T),grad); #local lsigma:=linalg[multiply](d,lgrad); #local sigma:=linalg[multiply](T,lsigma); #global w:=-linalg[multiply](sigma,grad)/2; U:=w*det_jac/2; map(eval,[convert(sigma,list),convert(grad,list),[U,w]])[] end proc: compute_fluxes:=proc(zelms::list) global nods,elems,e_velts,n_fluxes; local elms,dim,flx,flux,parms,e,velts,x,y,i,j,k,n,X,Y,Vx,Vy; if nargs=0 then elms:=[seq(i,i=1..nops(elems))] else elms:=zelms end if; flx:=element_fluxes([seq([X[i],Y[i]],i=1..3)],[Vx,Vy]); parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],Vx=velts[1],Vy=velts[2]]; k,dim:=0,nops(nods);n_fluxes:=[seq([0,0],i=1..dim)]; for e in elems do k:=k+1; if member(k,elms) then velts:=e_velts[k][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); flux:=subs(eval(parms)[],eval(flx)); n:=0; for i in e[1..3] do n:=n+1; n_fluxes[i]:=[seq(n_fluxes[i][j]+flux[n][j],j=1..2)] end do end if end do; n_fluxes:=simplify(n_fluxes,assume=positive); if nargs>1 then flx:=['q[x]','q[y]']: print(`Nodal fluxes (n_fluxes): `*flx);print(n_fluxes[]) end if; NULL end proc: element_fluxes:=proc(e_nods::list(list),e_vels::list) local i,q,b,c; q:=[];element_geometry(e_nods,b,c); for i to 3 do matrix([[b[i]/2,0],[0,c[i]/2]]);convert(e_vels,vector); q:=[q[],convert(linalg[multiply](%%,%),list)]; end do; eval(q) end proc: element_geometry:=proc(e_nods::list(list),b::evaln,c::evaln) local x,y,i,j,k,eb,ec; x,y:=seq([seq(e_nods[i][j],i=1..3)],j=1..2); for i to 3 do j,k:=loop_123(i); eb[i]:=y[j]-y[k]; ec[i]:=x[k]-x[j]; end do; if nargs>1 then b:=eval([seq(eb[i],i=1..3)]); c:=eval([seq(ec[i],i=1..3)]) end if; eval((eb[1]*ec[2]-eb[2]*ec[1])/2) end proc: loop_123:=proc(i::posint) local j,k; if i<3 then j:=i+1 else j:=1 end if; if j<3 then k:=j+1 else k:=1 end if; j,k end proc: end module: G_cgt_fem:=module() description "Graphics for 2D Potential Flow Analysis with CGT Finite Elements"; export plot_problem_data,plot_mesh,\ animate_3D_rot,animate_3D_def,\ animate_velocities,plot_velocities,\ plot_fluxes,plot_contour_lines,\ plot_line_potentials; local rotate_mesh,deform_mesh,rgb_color,\ stress_arrow,crossing_point,loop_123,\ mesh_line_intersection,e_side_cross_line,\ point_on_line; option package; plot_problem_data:=proc(llc::list,urc::list) global tcase,control,nods,elems,mat_props,bdr_conds,b_velts,p_srcs,e_srcs; local xlgth,ylgth,lgth,v,view,elem,n,g,x,y,xg,yg,i,j; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); Plotter:-start_2D_plot();Plotter:-draw_title(`Problem Data`); for elem in elems do; n:=elem[1..3];g:=elem[4]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; Plotter:-new_color(0,1,0);Plotter:-draw_text([xg,yg],convert(g,string),[1,1],10); Plotter:-new_color();Plotter:-move_pen([x[3],y[3]],0); seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3); end do; Plotter:-new_color(1,0,0); for n in bdr_conds do xg,yg:=nods[n[1]][1],nods[n[1]][2]; Plotter:-new_color();Plotter:-draw_symbol([xg,yg],3);Plotter:-new_color(1,0,0); Plotter:-draw_text([xg,yg],convert(evalf(n[2],3),string),[1,1],10); end do; Plotter:-new_color(1,0,1); if member("boundary velocities",control) then for n in b_velts do for i to 2 do xg,yg:=nods[n[i]][1],nods[n[i]][2]; Plotter:-new_color();Plotter:-draw_symbol([xg,yg],3);Plotter:-new_color(1,0,1); Plotter:-draw_text([xg,yg],convert(evalf(n[3],3),string),[2*i-3,2*i-3],10) end do end do end if; Plotter:-new_color(1,.7,0); if member("point sources",control) then for n in p_srcs do xg,yg:=nods[n[1]][1],nods[n[1]][2]; Plotter:-new_color();Plotter:-draw_symbol([xg,yg],3);Plotter:-new_color(1,.7,0); Plotter:-draw_text([xg,yg],convert(evalf(n[2],3),string),[1,-1],10); end do end if; if member("element sources",control) then for n in e_srcs do x,y:=seq([seq(nods[elems[n[1]][i]][k],i=1..3)],k=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; Plotter:-draw_text([xg,yg],convert(evalf(n[2],3),string),[-1,-1],10) end do end if; if nargs=2 then view:=args[1..2] else view:=NULL; [[1,0,0],[1,0,1],[1,.7,0],[0,1,0]]; Plotter:-draw_label(["Potential","Velocity","Source","Material"],%) end if; Plotter:-show_2D_plot(view); end proc: plot_mesh:=proc(llc::list,urc::list) global tcase,control,nods,elems,mat_props,bdr_conds,b_velts; local view,node,elem,n,x,y,i,j,xg,yg; Plotter:-start_2D_plot();Plotter:-draw_title(`Finite Element Mesh`); for elem in elems do; n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; member(elem,elems,'i');Plotter:-new_color(1,0,0); Plotter:-draw_text([xg,yg],convert(i,string),[0,0]);Plotter:-new_color(); x:=[seq(0.7*(x[i]-xg)+xg,i=1..3)];y:=[seq(0.7*(y[i]-yg)+yg,i=1..3)]; Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; Plotter:-new_color(0,0,1); for node in nods do member(node,nods,'i'); Plotter:-draw_text([node[1],node[2]],convert(i,string)) end do; if nargs=2 then view:=args[1..2] else view:=NULL; Plotter:-draw_label(["Node Number","Element Number"],[[0,0,1],[1,0,0]]) end if; Plotter:-show_2D_plot(view); end proc: animate_3D_rot:=proc(zv::list,llc::list,urc::list) global nods; local view,alpha,frms,xr,yr,i,j,zb; Plotter:-start_3D_animation(); Plotter:-draw_title(convert(readstat("plot title: "),string)); xr,yr:=seq(add(nods[i][j],i=1..nops(nods))/nops(nods),j=1..2); frms:=80;alpha:=2.0*Pi/frms;zb:=0; for i to frms do rotate_mesh(zv,[xr,yr],alpha*i,zb); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_3D_animation(view); end proc: rotate_mesh:=proc(zv::list,rotc::list,alpha::algebraic,zb::algebraic) global elems,nods; local nval,csa,sna,elem,n,h,x,y,i,j,k,xx,yy,xr,yr; if nops(zv)=nops(nods) then nval:=true else nval:=false end if; xr,yr:=rotc[]; csa,sna:=evalf(cos(alpha)),evalf(sin(alpha)); Plotter:-new_color(1,0,0);Plotter:-new_pen_style();Plotter:-new_thickness(1); Plotter:-new_plot_style();Plotter:-new_plot_shading(1);j:=0; for elem in elems do n:=elem[1..3];j:=j+1; if nval then h:=[seq(zv[n[i]],i=1..3)] else h:=[seq(zv[j],i=1..3)] end if; x,y:=seq([seq(nods[n[i]][k],i=1..3)],k=1..2); if nargs>1 then xx:=[seq(csa*(x[i]-xr)+sna*(y[i]-yr),i=1..3)]; yy:=[seq(-sna*(x[i]-xr)+csa*(y[i]-yr),i=1..3)]; x,y:=xx,yy end if; Plotter:-new_color(1,0,0);Plotter:-plot_element(x,y,h); Plotter:-new_color(1,1,0);Plotter:-plot_element(x,y,[zb,zb,zb]) end do; NULL end proc: animate_3D_def:=proc(zv::list,llc::list,urc::list) local view,frms,i,j; Plotter:-start_3D_animation(); Plotter:-draw_title(convert(readstat("plot title: "),string)); frms:=80; for i to frms do j:=2.0*i/frms;if j>1 then j:=2-j end if; deform_mesh([seq(zv[k]*j,k=1..nops(zv))]); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_3D_animation(view); end proc: deform_mesh:=proc(zv::list,rgb::list) global elems,nods; local nval,elem,n,h,x,y,i,j,k; if nops(zv)=nops(nods) then nval:=true else nval:=false fi; j:=0;Plotter:-new_pen_style();Plotter:-new_plot_style();Plotter:-new_plot_shading(); for elem in elems do n:=elem[1..3];j:=j+1; if nval then h:=[seq(zv[n[i]],i=1..3)] else h:=[seq(zv[j],i=1..3)] end if; if nargs=1 then x,y:=seq([seq(nods[n[i]][k],i=1..3)],k=1..2); Plotter:-new_color();Plotter:-new_thickness(1); Plotter:-plot_element(x,y,h);Plotter:-plot_element(x,y,[0,0,0]) else x:=[seq(nods[n[i]][1]+h[i][1],i=1..3)]; y:=[seq(nods[n[i]][2]+h[i][2],i=1..3)]; Plotter:-new_color(rgb[]);Plotter:-new_thickness();Plotter:-plot_element(x,y) end if; end do; NULL end proc: plot_contour_lines:=proc(nlevels::list,lvls::list,llc::list,urc::list) global elems,nods; local elem,level,levels,line,view,x,y,p_in,p_out,d,dm,n,z,i,j,k,cols,a,d0; Plotter:-start_2D_plot(); Plotter:-draw_title(convert(readstat("plot title: "),string)); if nargs>1 then levels:=lvls else d:=evalf((max(nlevels[])-min(nlevels[]))/9); levels:=[seq(min(nlevels[])+i*d,i=0..9)] end if; Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for elem in elems do; n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; Plotter:-new_pen_style();Plotter:-new_thickness(3); d0:=min(levels[]);dm:=max(levels[])-d0;cols:=[]; if abs(dm)<10^(5-Digits) then RETURN(`identical levels`) end if; for level in levels do d:=(level-d0)/dm;a:=evalf(cos(d*Pi)); Plotter:-new_color(rgb_color(a));cols:=[cols[],[rgb_color(a)]]; line:=[]; for elem in elems do n:=elem[1..3];z:=seq(nlevels[n[i]],i=1..3); for i to 3 do j,k:=loop_123(i); if z[i]<=level and z[j]>=level then p_in:=crossing_point([n[j],n[i]],[z[j],z[i]],level); if z[k][] then for i to nops(line) do Plotter:-move_pen(line[i][1],0);Plotter:-move_pen(line[i][2],1) end do end if end do; if nargs>2 then view:=args[3..4] else view:=NULL; [seq(convert(evalf(levels[i],2),string),i=1..nops(levels))],cols; Plotter:-draw_label(%) end if; Plotter:-show_2D_plot(view); end proc: crossing_point:=proc(n::list,z::list,level::algebraic) global nods; local t,i; if abs(z[1]-z[2])<10^(1-Digits) then t:=[1] else t:=[evalf((level-z[2])/(z[1]-z[2]))] end if; t:=[t[],evalf(1-t[1])]; [add(t[i]*nods[n[i]][1],i=1..2),add(t[i]*nods[n[i]][2],i=1..2)] end proc: loop_123:=proc(i::posint) local j,k; if i<3 then j:=i+1 else j:=1 end if; if j<3 then k:=j+1 else k:=1 end if; j,k end proc: rgb_color:=proc(a) local x,n,m,r,g,b; x:=min(a,1);x:=max(x,-1); n:=1.3;m:=1.5; r:=1/(1+n**2*(x+1)**2); g:=1/(1+m**2*x**2); b:=1/(1+n**2*(x-1)**2); r,g,b end proc: plot_fluxes:=proc(elms::list,llc::list,urc::list) global nods,elems,n_velts,n_fluxes; local e,x,y,mxv,scl,n,v,i,j,lab,view,xlgth,ylgth,lgth; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); if nargs>0 then Cgt_fem:-compute_fluxes(elms) else Cgt_fem:-compute_fluxes([seq(i,i=1..nops(elems))]) end if; Plotter:-start_2D_plot();Plotter:-draw_title(`Nodal Flux Through Element Sides`); Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for e in elems do n:=e[1..3];x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; mxv:=max(seq(n_fluxes[i][],i=1..nops(n_fluxes)));scl:=lgth/mxv/2; for i to nops(nods) do if evalf(abs(add(n_fluxes[i][j],j=1..2)))>10^(5-Digits) then n:=nods[i];v:=evalf([n_fluxes[i][1]*scl,n_fluxes[i][2]*scl]); if v[1]>0 then Plotter:-new_color(1,0,1) else Plotter:-new_color(0,1,0) end if; Plotter:-new_arrow(n,[v[1],0],v[1]*0.25,0,0); Plotter:-new_arrow(n,[-v[1],0],v[1]*0.25,0,0); if v[2]>0 then Plotter:-new_color(1,0,1) else Plotter:-new_color(0,1,0) end if; Plotter:-new_arrow(n,[0,v[2]],v[2]*0.25,0,0); Plotter:-new_arrow(n,[0,-v[2]],v[2]*0.25,0,0) end if end do; lab:=cat("Max Value: ",convert(evalf(mxv,2),string)); if nargs>1 then view:=args[2..3] else view:=NULL; Plotter:-draw_label(["Inflow",lab,"Outflow"],[[0,1,0],[0,0,0],[1,0,1]]) end if; Plotter:-show_2D_plot(view); end proc: animate_velocities:=proc(vscale,llc::list,urc::list) global nods,elems,n_velts; local view,frms,e,i,j,k,l,x,y,mxv,scl,scv,n,xlgth,ylgth,lgth; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); mxv:=max(seq(n_velts[i][1]^2+n_velts[i][2]^2,i=1..nops(n_velts))); scv:=1;if nargs>0 then scv:=max(scv,args[1]) end if; scl:=evalf(lgth/sqrt(mxv))/scv; frms:=50;Plotter:-start_2D_animation(); Plotter:-draw_title(`Animated Flow Velocities`); for i to frms do j:=2.0*i/frms;if j>1 then j:=2-j end if; Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for e in elems do n:=e[1..3]; x,y:=seq([seq(nods[n[k]][l],k=1..3)],l=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[k],y[k]],1),k=1..3) end do; evalf([seq([seq(n_velts[k][l]*scl*j,l=1..2)],k=1..nops(nods))]),[.4,.9,.6]; deform_mesh(%); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_2D_animation(view); end proc: plot_velocities:=proc(vscale,llc::list,urc::list) global nods,elems,n_velts; local elem,x,y,mxv,mnv,dv,scl,scv,n,v,i,j,lab,view,xlgth,ylgth,lgth,th,a,xg,yg; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); Plotter:-start_2D_plot();Plotter:-draw_title(`Flow Velocity Field`); Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for elem in elems do n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; scv:=1;if nargs>0 then scv:=max(scv,args[1]) end if; mxv:=max(seq(n_velts[i][1]^2+n_velts[i][2]^2,i=1..nops(n_velts))); mxv:=sqrt(mxv);scl:=evalf(lgth/mxv)/scv; mnv:=min(seq(n_velts[i][1]^2+n_velts[i][2]^2,i=1..nops(n_velts))); mnv:=sqrt(mnv)*scl;dv:=evalf(abs(lgth-mnv))/scv; for i to nops(elems) do n:=elems[i][1..3];x,y:=seq([seq(nods[n[k]][j],k=1..3)],j=1..2); xg:=add(x[k],k=1..3)/3;yg:=add(y[k],k=1..3)/3; n:=[xg,yg];v:=evalf([e_velts[i][1]*scl,e_velts[i][2]*scl]); th:=evalf(sqrt(v[1]^2+v[2]^2)); if dv<10^(1-Digits) then a:=-1 else a:=evalf(cos(abs(th-mnv)/dv*Pi)) end if; Plotter:-new_color(rgb_color(a));Plotter:-new_arrow(n,v,th*0.25,th*0.5,1/3) end do; for i to nops(nods) do n:=nods[i];v:=evalf([n_velts[i][1]*scl,n_velts[i][2]*scl]); th:=evalf(sqrt(v[1]^2+v[2]^2)); if dv<10^(1-Digits) then a:=-1 else a:=evalf(cos(abs(th-mnv)/dv*Pi)) end if; Plotter:-new_color(rgb_color(a));Plotter:-new_arrow(n,v,th*0.25,th*0.5,1/3) end do; lab:=cat("Maximum Value: ",convert(evalf(mxv,2),string)); if nargs>1 then view:=args[2..3] else view:=NULL;Plotter:-draw_label([lab],[[1,0,0]]) end if; Plotter:-show_2D_plot(view); end proc: plot_line_potentials:=proc(zval::list) global n_potls; local i,zv,curve,curves,line,lines,labels,labcolors,going,rep; curves,lines,labels,labcolors:=[],[],[],[];going:=true; while going do rep:=evalf(readstat(" enter line end-points as X1,Y1,X2,Y2 or 0 to finish: ")); if rep=0 then going:=false else lines:=[lines[],[[rep[1..2]],[rep[3..4]]]]; rep:=convert(readstat(" line label: "),string); labels:=[labels[],rep]; rep:=readstat(" color as R,G,B: "); labcolors:=[labcolors[],[rep]] end if end do; rep:=convert(readstat(" plot title: "),string); Plotter:-start_2D_plot();Plotter:-new_axes(1);Plotter:-draw_title(rep); if nargs=0 then zv:=n_potls else zv:=zval end if; Plotter:-new_thickness(2); for line in lines do curves:=[curves[],mesh_line_intersection(line,zv)] end do; for curve in curves do member(curve,curves,'i'); if nops(curve)>1 then Plotter:-new_color(labcolors[i][]); seq(Plotter:-draw_symbol(curve[i],1),i=1..nops(curve)); Plotter:-move_pen(curve[1],0);seq(Plotter:-move_pen(curve[i],1),i=2..nops(curve)) else labels:=subsop(i=NULL,labels);labcolors:=subsop(i=NULL,labcolors) end if end do; Plotter:-draw_label(labels,labcolors); Plotter:-show_2D_plot(); end proc: mesh_line_intersection:=proc(line::listlist,zval::list) global nods,elems; local element,i,j,k,curve,ni,nj,l_line,x,y,cross,l_xy,n_i,n_j,h,first; curve:=[];ni:=[line[1],line[2],0];nj:=[line[1],line[2],0]; l_line:=sqrt((line[2][1]-line[1][1])^2+(line[2][2]-line[1][2])^2); for element in elems do for i to 3 do j:=i+1;if j>3 then j:=1 end if; ni[3]:=nods[element[i]];nj[3]:=nods[element[j]]; cross,x,y,n_j:=e_side_cross_line(ni,nj); if cross then l_xy:=sqrt((x-line[1][1])^2+(y-line[1][2])^2); if point_on_line(x,y,line,l_xy/l_line) then if curve=[] then first:=l_xy else first:=min(first,l_xy) end if; n_i:=1-n_j;h:=zval[element[i]]*n_i+zval[element[j]]*n_j; if not member([l_xy,h],curve) then k:=nops(curve);curve:=[curve[],[]]; while k>0 and l_xy=0 and ratio<=1 and abs(x-line[1][1]-ratio*(line[2][1]-line[1][1]))<10^(1-Digits) and abs(y-line[1][2]-ratio*(line[2][2]-line[1][2]))<10^(1-Digits) then true else false end if end proc: end module: Cst_fem:=module() description "2D Linear Elastic Analysis with CST Finite Elements"; export cst_fem,read_save_data; local principal_stresses,shape_f,lder_shape_f,jacob,gder_shape_f,\ initialize_plane_condition,initialize_global_matrices,\ assemble_global_matrices,force_exact_boundary_conditions,\ compute_displacements,compute_reactions,compute_stresses,\ element_k,set_up_d_b,assemble_stiffness,assemble_point_forces,\ assemble_self_weight,element_p,skew_constraints,rotation_matrix,\ constrained_eqtn,element_stresses,read_data,read_line,\ line_contents,save_data,rotate_mesh,deform_mesh,rgb_color; option package; cst_fem:=proc() # *** main program *** global nods,g_k,g_p,disp; local cpu_time,talk,u; cpu_time:=time(); if nargs=0 then talk:=NULL else talk:=true end if; initialize_plane_condition(talk); initialize_global_matrices(talk); assemble_global_matrices(talk); force_exact_boundary_conditions(talk); u:=linalg[linsolve](g_k,g_p); compute_displacements(u,talk); compute_reactions(u,talk); compute_stresses(talk); cpu_time:=convert(time()-cpu_time,string); print(`Cpu time: `||cpu_time||` seconds`); NULL end proc: initialize_plane_condition:=proc() global plane_case,mat_props; local mat,i; if plane_case="plane strain" then i:=0; for mat in mat_props do mat[2]:=mat[2]/(1+mat[2]); mat[1]:=mat[1]*(1-mat[2]**2); i:=i+1;mat_props[i]:=mat end do; else plane_case:="plane stress"; end if; if nargs<>0 then print(`Plane`*convert(plane_case[7..12],name)*`analysis`) end if; NULL end proc: initialize_global_matrices:=proc() global nods,g_k,g_p; local dim,i,j,p; dim:=2*nops(nods); g_k:=array(symmetric,1..dim,1..dim,[seq([seq(0,j=1..i)],i=1..dim)]); g_p:=array([seq(0,j=1..dim)]); if nargs<>0 then print(`Initializing global matrices:`);print(g_k,g_p) end if; NULL end proc: assemble_global_matrices:=proc() global control; local c; assemble_stiffness(args); for c in control do if c="point forces" then assemble_point_forces(args) elif c="self weight" then assemble_self_weight(args) end if end do; NULL end proc: force_exact_boundary_conditions:=proc() global reacts,bdr_conds,nods,g_k,g_p; local dim,cond,eqtn,n,i,j; dim:=2*nops(nods); skew_constraints(); for cond in bdr_conds do eqtn:=constrained_eqtn(cond); n:=2*(cond[1]-1)+max(1,eqtn); for i to dim do g_p[i]:=g_p[i]-cond[3]*g_k[i,n] end do; for i to dim do g_k[n,i]:=0 end do; g_k[n,n]:=1; g_p[n]:=cond[3] end do; if nargs<>0 then print(`Forcing exact boundary conditions:`); print(simplify(g_k,assume=positive)); print(simplify(g_p,assume=positive)) end if; NULL end proc: compute_displacements:=proc(u::{list,array}) global disp,bdr_conds; local cond,TT,n,i; disp:=[seq([u[2*i-1],u[2*i]],i=1..nops(nods))]; for cond in bdr_conds do if constrained_eqtn(cond)=0 then TT:=rotation_matrix(cond[2]*Pi/180); n:=2*cond[1]; disp[cond[1]]:=linalg[multiply](TT,array([u[n-1],u[n]])) end if end do; disp:=simplify(disp); if nargs>1 then print(`Nodal displacements (disp):`*['u[x]','u[y]']);print(disp[]) end if; NULL end proc: compute_reactions:=proc(u::{list,array}) global reacts,nods,g_p; local i,j,dim,c,s,react,eqtn,fx,fy,r,p; i,fx,fy,dim:=0,0,0,2*nops(nods); for react in reacts do i:=i+1;r:=add(react[3][j]*u[j],j=1..dim); eqtn:=constrained_eqtn(react); if eqtn=1 or eqtn=2 then reacts[i]:=[seq(react[j],j=1..2),[r]]; if eqtn=1 then fx:=fx+r else fy:=fy+r fi; else react[2]*Pi/180;s:=sin(%);c:=cos(%%); reacts[i]:=[react[1],react[2],[c*r,s*r]]; fx,fy:=fx+c*r,fy+s*r; end if end do; for p in forces do fx:=fx+p[2];fy:=fy+p[3] end do; reacts:=simplify(reacts,assume=positive); fx:=simplify(fx,assume=positive); fy:=simplify(fy,assume=positive); if nargs>1 then print(`Reactions (reacts):`*['node','direction',['reaction']]); print(reacts[]);fx,fy:=convert(fx,string),convert(fy,string); print(`Equilibrium: `*'Sigma*F[x]'*`= `||fx,` and `*'Sigma*F[y]'*`= `||fy) end if; NULL end proc: compute_stresses:=proc() global mat_props,nods,disp,e_sigma,n_sigma,total_strain_energy,symblc_prob; local dim,sig,sigma,parms,hit,e,mat,x,y,v,i,j,X,Y,V,E,nu,th; sig:=element_stresses([seq([X[i],Y[i]],i=1..3)],[E,nu,th],[seq(V[i],i=1..6)]); parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],E=mat[1],nu=mat[2],th=mat[4]]; parms:=[parms[],seq(V[i]=v[i],i=1..6)]; dim:=nops(nods);e_sigma:=[];n_sigma:=[seq([0,0,0,0],i=1..dim)]; hit:=[seq(0,i=1..dim)];total_strain_energy:=0; for e in elems do mat:=mat_props[e[4]][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); v:=[seq(seq(disp[e[i]][j],j=1..2),i=1..3)]; sigma:=subs(eval(parms)[],eval(sig)); total_strain_energy:=total_strain_energy+sigma[3]; e_sigma:=[e_sigma[],sigma[1..3]]; for i in e[1..3] do [seq(n_sigma[i][j]+sigma[1][j],j=1..3),n_sigma[i][4]+sigma[4]]; n_sigma[i]:=%; hit[i]:=hit[i]+1 end do end do; n_sigma:=[seq([seq(n_sigma[i][j]/hit[i],j=1..4)],i=1..dim)]; for i to dim do n_sigma[i]:=[principal_stresses([seq(n_sigma[i][j],j=1..3)])[],n_sigma[i][4]]; end do; e_sigma:=simplify(e_sigma); n_sigma:=simplify(n_sigma,assume=positive); total_strain_energy:=simplify(total_strain_energy,assume=positive); if nargs<>0 then sig:=[['sigma[x]','sigma[y]','sigma[xy]']]: sig:=[sig[],['sigma[I]','sigma[II]','angle'],`strain energy`]: print(`Element stresses (e_sigma): `*sig);print(e_sigma[]); sig:=['sigma[I]','sigma[II]','angle',`strain energy density`]: print(`Nodal stresses (n_sigma): `*sig);print(n_sigma[]); print(`Total strain energy: `);print(total_strain_energy) end if; NULL end proc: shape_f:=proc(nat_crds::list) local s,t,N; option remember; N:=[1-s-t,s,t]; # nodes at: 1->(0,0) 2->(1,0) 3->(0,1) if nargs<>0 then N:=subs(s=nat_crds[1],t=nat_crds[2],N) end if; eval(N) end proc: lder_shape_f:=proc(nat_crds::list) local N,s,t,lder,i; option remember; N:=shape_f([s,t]); lder:=[seq([diff(N[i],s),diff(N[i],t)],i=1..3)]; if nargs<>0 then lder:=subs(s=nat_crds[1],t=nat_crds[2],lder) end if; eval(lder) end proc: jacob:=proc(e_nods::list(list),nat_crds::list) local nod_crds,N,s,t,x,y,jac,i; option remember; nod_crds:=map(i->i[],e_nods); N:=shape_f([s,t]); N:=seq(array([[N[i],0],[0,N[i]]]),i=1..3); N:=linalg[augment](N[1],N[2],N[3]);linalg[multiply](N,nod_crds); x,y:=%[1],%[2]; jac:=linalg[transpose](linalg[jacobian]([x,y],[s,t])); if nargs>1 then jac:=subs(s=nat_crds[1],t=nat_crds[2],eval(jac)) end if; eval(jac); end proc: gder_shape_f:=proc(e_nods::list(list),det_jac::evaln,nat_crds::list) local jac,i_jac,s,t,gder,i; option remember; jac:=jacob(e_nods,[s,t]); det_jac:=linalg[det](jac); i_jac:=linalg[inverse](jac); gder:=[seq(linalg[multiply](i_jac,lder_shape_f([s,t])[i]),i=1..3)]; if nargs>2 then gder:=subs(s=nat_crds[1],t=nat_crds[2],eval(gder)) end if; eval(gder) end proc: element_k:=proc(e_nods::list(list),e_props::list) local s,t,d,b,det_jac,kernel; d,b,det_jac:=set_up_d_b(e_nods,e_props,[s,t]); kernel:=linalg[multiply](linalg[transpose](b),d,b)*det_jac; evalm(int(int(kernel,t=0..1-s),s=0..1)*e_props[3]) end proc: set_up_d_b:=proc(e_nods::list(list),e_props::list,nat_crds::list) local E,nu,d,gder,s,t,b,det_jac,i; option remember; E,nu:=e_props[1],e_props[2]; d:=array([[1,nu,0],[nu,1,0],[0,0,(1-nu)/2]])*E/(1-nu**2); gder:=gder_shape_f(e_nods,det_jac,[s,t]); b:=seq(array([[gder[i][1],0],[0,gder[i][2]],[gder[i][2],gder[i][1]]]),i=1..3); b:=linalg[augment](b[1],b[2],b[3]); if nargs>2 then det_jac:=subs(s=nat_crds[1],t=nat_crds[2],det_jac); b:=subs(s=nat_crds[1],t=nat_crds[2],eval(b)) end if; map(eval,[d,b,det_jac])[] end proc: assemble_stiffness:=proc() global nods,elems,mat_props,g_k; local e_k,ek,parms,e,n,mat,i,j,x,y,X,Y,E,nu,th; e_k:=element_k([seq([X[i],Y[i]],i=1..3)],[E,nu,th]); parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],E=mat[1],nu=mat[2],th=mat[4]]; for e in elems do mat:=mat_props[e[4]][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); ek:=simplify(subs(eval(parms)[],eval(e_k)),assume=positive); n:=map(proc(i) 2*i-1,2*i end, e[1..3]); for i to 6 do for j to i do g_k[n[i],n[j]]:=g_k[n[i],n[j]]+ek[i,j] end do end do; if nargs<>0 then print('`Assembling element`'*e*'`:`'); print(simplify(g_k,assume=positive)) end if end do; NULL end proc: assemble_point_forces:=proc() global forces,g_p; local p,i; for p in forces do i:=2*p[1]-1; g_p[i]:=g_p[i]+p[2]; g_p[i+1]:=g_p[i+1]+p[3] end do; if nargs<>0 then print(`Assembling point forces:`);print(g_p) end if; NULL end proc: assemble_self_weight:=proc() global nods,elems,mat_props,g_p; local e_p,ep,parms,e,n,mat,i,j,x,y,X,Y,fy; e_p:=element_p([seq([X[i],Y[i]],i=1..3)],[fy,th]); parms:=[seq(X[i]=x[i],i=1..3),seq(Y[i]=y[i],i=1..3)]; parms:=[parms[],fy=mat[3],th=mat[4]]; for e in elems do mat:=mat_props[e[4]][]; x,y:=seq([seq(nods[e[i]][j],i=1..3)],j=1..2); ep:=simplify(subs(eval(parms)[],eval(e_p)),assume=positive); n:=map(proc(i) 2*i-1,2*i end, e[1..3]); for i to 6 do g_p[n[i]]:=g_p[n[i]]+ep[i] end do end do; if nargs<>0 then print(`Assembling self weight:`); print(simplify(g_p,assume=positive)) end if; NULL end proc: element_p:=proc(e_nods::list(list),e_props::list) local fy,th,bf,det_jac,s,t,N,i,kernel; fy,th:=e_props[1],e_props[2];bf:=array([0,fy]); det_jac:=linalg[det](jacob(e_nods,[s,t])); N:=shape_f([s,t]); N:=seq(array([[N[i],0],[0,N[i]]]),i=1..3); N:=linalg[augment](N[1],N[2],N[3]); kernel:=convert(linalg[multiply](linalg[transpose](N),bf),list); map(k->int(int(k*det_jac,t=0..1-s),s=0..1)*th,kernel,det_jac,s,t,th) end proc: skew_constraints:=proc() ### only one constraint/node is possible global reacts,bdr_conds,g_k,g_p; local dim,cond,eqtn,alpha,TT,n,cols,diag,j; reacts:=[];dim:=linalg[coldim](g_k); for cond in bdr_conds do eqtn:=constrained_eqtn(cond); if eqtn=0 then TT:=rotation_matrix(cond[2]*Pi/180); n:=2*cond[1]; linalg[submatrix](g_k,1..linalg[coldim](g_k),n-1..n); cols:=linalg[multiply](%,TT); linalg[submatrix](cols,n-1..n,1..2); diag:=linalg[multiply](linalg[transpose](TT),%); linalg[copyinto](diag,cols,n-1,1); linalg[copyinto](cols,g_k,1,n-1); # note g_k is declared as symmetric! array([seq(g_p[j],j=n-1..n)]); linalg[multiply](linalg[transpose](TT),%); g_p[n-1],g_p[n]:=%[1],%[2] end if end do; for cond in bdr_conds do eqtn:=constrained_eqtn(cond); n:=2*(cond[1]-1)+max(1,eqtn); reacts:=[reacts[],[cond[1],cond[2],[seq(g_k[n,j],j=1..dim)]]] end do; NULL end proc: rotation_matrix:=proc(alpha::algebraic,s::evaln,c::evaln) local sn,cn; sn,cn:=sin(alpha),cos(alpha); if nargs>1 then s,c:=sn,cn fi; array([[cn,-sn],[sn,cn]]) end proc: constrained_eqtn:=proc(cond::list) local eqtn,dir; eqtn:=0; dir:=convert(cond[2],string); if dir="x" then eqtn:=1 elif dir="y" then eqtn:=2 end if; eqtn end proc: element_stresses:=proc(e_nods::list(list),e_props::list,u::list) local d,b,det_jac,sigma,p_sigma,E,nu,th,w,U; d,b,det_jac:=set_up_d_b(e_nods,e_props); sigma:=linalg[multiply](d,b,u); p_sigma:=principal_stresses(sigma); E,nu,th:=e_props[1],e_props[2],e_props[3]; w:=(p_sigma[1]**2+p_sigma[2]**2-2*nu*p_sigma[1]*p_sigma[2])/(2*E); U:=w*det_jac/2*th; map(eval,[sigma,p_sigma,U,w]) end proc: principal_stresses:=proc(sigma::{list,array}) local c,sig,r; c:=(sigma[1]+sigma[2])/2; sig:=(sigma[1]-sigma[2])/2; r:=sqrt(sig**2+sigma[3]**2); map(eval,[c+r,c-r,arctan(sigma[3],sig)*90/Pi]) end proc: read_save_data:=proc() local reply; reply:=0; while reply<>y and reply<>n do reply:=readstat("read data from a file (y/n) ? ") end do; if reply=y then read_data() end if; reply:=0; while reply<>y and reply<>n do reply:=readstat("save data into a file (y/n) ? ") end do; if reply=y then save_data() end if; NULL end proc: read_data:=proc() global tcase,plane_case,control,nods,elems,forces,mat_props,bdr_conds; local reply,fle,nod,elem,line,lne,i; reply:=readstat("file name: "); fle:=fopen(reply,READ,TEXT); control,nods,elems,forces,mat_props,bdr_conds:=seq([],i=1..6); line:=sscanf(readline(fle),"%s"); while line<>["*end*"] do if line=["*control*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%s"); if substring(lne,1..5)="title" then tcase:=lne[6..length(lne)] elif substring(lne,1..5)="plane" then plane_case:=substring(lne,1..12) elif substring(lne,1..12)="point forces" then control:=[control[],"point forces"] elif substring(lne,1..11)="self weight" then control:=[control[],"self weight"] fi end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*nodes*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a"); if line<>0 and line<>[] then line_contents(line,[nods[],[seq(line[i],i=2..3)]],nods) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*elements*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a%a%a"); if line<>0 and line<>[] then line_contents(line,[elems[],[seq(line[i],i=2..5)]],elems) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*forces*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a"); if line<>0 and line<>[] then line_contents(line,[forces[],[seq(line[i],i=1..3)]],forces) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*constraints*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a"); if line<>0 and line<>[] then line_contents(line,[bdr_conds[],[seq(line[i],i=1..3)]],bdr_conds) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if elif line=["*materials*"] then print(line[]); while line<>0 and line<>[] do lne,line:=read_line(fle,"%a%a%a%a%a"); if line<>0 and line<>[] then line_contents(line,[mat_props[],[seq(line[i],i=2..5)]],mat_props) end if end do; if line<>0 and line<>[] then line:=sscanf(lne,"%s") end if else line:=sscanf(readline(fle),"%s") end if end do; fclose(fle); NULL end proc: read_line:=proc(fle::integer,format::string) local lne; lne:=readline(fle); lne,sscanf(lne,format) end proc: line_contents:=proc(line::{list,integer},lst::list,var::evaln) if line<>0 and line<>[] then var:=lst else NULL end if end proc: save_data:=proc() global tcase,plane_case,control,nods,elems,forces,mat_props,bdr_conds; local reply,file,i; reply:=readstat("file name: "); file:=fopen(reply,WRITE,TEXT); fprintf(file,"\n*control* [title,plane stress/strain,"); fprintf(file,"point forces,self weight]\n"); fprintf(file,"%s%s\n%s\n","title",tcase,plane_case); map(proc(i,file) fprintf(file,"%s\n",i) end,control,file); fprintf(file,"\n*materials* [material,Young modulus,Poisson coef,"); fprintf(file,"specific weight,thickness]\n"); seq(fprintf(file,"%a %a %a %a %a\n",i,mat_props[i][]),i=1..nops(mat_props)); fprintf(file,"\n*nodes* [node,x,y]\n"); seq(fprintf(file,"%a %a %a\n",i,nods[i][]),i=1..nops(nods)); fprintf(file,"\n*elements* [element,node1,node2,node3,material]\n"); seq(fprintf(file,"%a %a %a %a %a\n",i,elems[i][]),i=1..nops(elems)); fprintf(file,"\n*constraints* [node,direction"); fprintf(file,"(x,y or angle measured from x),displacement]\n"); map(proc(i,file) fprintf(file,"%a %a %a\n",i[]) end,bdr_conds,file); if nops(forces)>0 then fprintf(file,"\n*forces* [node,fx,fy]\n"); map(proc(i,file) fprintf(file,"%a %a %a\n",i[]) end,forces,file) end if; fprintf(file,"\n*end*\n"); fclose(file); NULL end proc: end module: G_cst_fem:=module() description "Graphics for 2D Linear Elastic Analysis with CST Finite Elements"; export plot_problem_data,plot_mesh,\ animate_3D_rot,animate_3D_def,\ animate_displacements,plot_displacements,\ plot_stresses,plot_contour_lines; local rotate_mesh,deform_mesh,rgb_color,\ stress_arrow,crossing_point,loop_123; option package; plot_problem_data:=proc(llc::list,urc::list) global tcase,plane_case,control,nods,elems,forces,mat_props,bdr_conds; local xlgth,ylgth,lgth,mxv,scl,v,th,view,elem,n,g,x,y,xg,yg,i,j; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); Plotter:-start_2D_plot();Plotter:-draw_title(`Problem Data`); for elem in elems do; n:=elem[1..3];g:=elem[4]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; Plotter:-new_color(0,1,0);Plotter:-draw_text([xg,yg],convert(g,string),[0,-1],10); Plotter:-new_color();Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3); end do; Plotter:-new_color(1,0,0); mxv:=max(seq(forces[i][2]^2+forces[i][3]^2,i=1..nops(forces))); scl:=evalf(lgth/sqrt(mxv)); for n in forces do v:=[n[2]*scl,n[3]*scl];th:=sqrt(v[1]^2+v[2]^2); Plotter:-new_arrow([nods[n[1]][1],nods[n[1]][2]],v,th*0.15,th*0.5,1/3) end do; Plotter:-new_color(1,0,1);lgth:=lgth*0.3; for n in bdr_conds do if convert(n[2],string)="x" or n[2]=0 then v:=[lgth,0] elif convert(n[2],string)="y" or n[2]=90 then v:=[0,lgth] else n[2]*evalf(Pi)/180;v:=[lgth*cos(%),lgth*sin(%)] end if; lgth*0.8;Plotter:-new_arrow([nods[n[1]][1],nods[n[1]][2]],v,%,%,1) end do; if nargs=2 then view:=args[1..2] else view:=NULL; Plotter:-draw_label(["Load","Constraint","Material"],[[1,0,0],[1,0,1],[0,1,0]]) end if; Plotter:-show_2D_plot(view); end proc: plot_mesh:=proc(llc::list,urc::list) global tcase,plane_case,control,nods,elems,forces,mat_props,bdr_conds; local view,node,elem,n,x,y,i,j,xg,yg; Plotter:-start_2D_plot();Plotter:-draw_title(`Finite Element Mesh`); for elem in elems do; n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; member(elem,elems,'i');Plotter:-new_color(1,0,0); Plotter:-draw_text([xg,yg],convert(i,string),[0,0]);Plotter:-new_color(); x:=[seq(0.7*(x[i]-xg)+xg,i=1..3)];y:=[seq(0.7*(y[i]-yg)+yg,i=1..3)]; Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; Plotter:-new_color(0,0,1); for node in nods do member(node,nods,'i'); Plotter:-draw_text([node[1],node[2]],convert(i,string)) end do; if nargs=2 then view:=args[1..2] else view:=NULL; Plotter:-draw_label(["Node Number","Element Number"],[[0,0,1],[1,0,0]]) end if; Plotter:-show_2D_plot(view); end proc: animate_3D_rot:=proc(zv::list,llc::list,urc::list) global nods; local view,alpha,frms,xr,yr,i,j; Plotter:-start_3D_animation(); Plotter:-draw_title(convert(readstat("plot title: "),string)); xr,yr:=seq(add(nods[i][j],i=1..nops(nods))/nops(nods),j=1..2); frms:=80;alpha:=2.0*Pi/frms; for i to frms do rotate_mesh(zv,[xr,yr],alpha*i); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_3D_animation(view); end proc: rotate_mesh:=proc(zv::list,rotc::list,alpha::algebraic) global elems,nods; local nval,csa,sna,elem,n,h,x,y,i,j,k,xx,yy,xr,yr,a,b,c; if nops(zv)=nops(nods) then nval:=true else nval:=false end if; if nargs>1 then xr,yr:=rotc[]; csa,sna:=evalf(cos(alpha)),evalf(sin(alpha)) end if; j:=0;Plotter:-new_pen_style();Plotter:-new_plot_style();Plotter:-new_plot_shading(); for elem in elems do n:=elem[1..3];j:=j+1; if nval then h:=[seq(zv[n[i]],i=1..3)] else h:=[seq(zv[j],i=1..3)] end if; x,y:=seq([seq(nods[n[i]][k],i=1..3)],k=1..2); if nargs>1 then xx:=[seq(csa*(x[i]-xr)+sna*(y[i]-yr),i=1..3)]; yy:=[seq(-sna*(x[i]-xr)+csa*(y[i]-yr),i=1..3)]; x,y:=xx,yy end if; Plotter:-new_color();Plotter:-new_thickness(1); Plotter:-plot_element(x,y,h);Plotter:-plot_element(x,y,[0,0,0]); a:=[x[1],x[1],x[2],x[2],x[2],x[2],x[3],x[3],x[3],x[3],x[1],x[1]]; b:=[y[1],y[1],y[2],y[2],y[2],y[2],y[3],y[3],y[3],y[3],y[1],y[1]]; c:=[h[1],0,0,h[2],h[2],0,0,h[3],h[3],0,0,h[1]]; Plotter:-new_color(1,0.8,0);Plotter:-new_thickness();Plotter:-plot_element(a,b,c) end do; NULL end proc: animate_3D_def:=proc(zv::list,llc::list,urc::list) local view,frms,i,j; Plotter:-start_3D_animation(); Plotter:-draw_title(convert(readstat("plot title: "),string)); frms:=80; for i to frms do j:=2.0*i/frms;if j>1 then j:=2-j end if; deform_mesh([seq(zv[k]*j,k=1..nops(zv))]); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_3D_animation(view); end proc: deform_mesh:=proc(zv::list,rgb::list) global elems,nods; local nval,elem,n,h,x,y,i,j,k,a,b,c; if nops(zv)=nops(nods) then nval:=true else nval:=false fi; j:=0;Plotter:-new_pen_style();Plotter:-new_plot_style();Plotter:-new_plot_shading(); for elem in elems do n:=elem[1..3];j:=j+1; if nval then h:=[seq(zv[n[i]],i=1..3)] else h:=[seq(zv[j],i=1..3)] end if; if nargs=1 then x,y:=seq([seq(nods[n[i]][k],i=1..3)],k=1..2); Plotter:-new_color();Plotter:-new_thickness(1); Plotter:-plot_element(x,y,h);Plotter:-plot_element(x,y,[0,0,0]); a:=[x[1],x[1],x[2],x[2],x[2],x[2],x[3],x[3],x[3],x[3],x[1],x[1]]; b:=[y[1],y[1],y[2],y[2],y[2],y[2],y[3],y[3],y[3],y[3],y[1],y[1]]; c:=[h[1],0,0,h[2],h[2],0,0,h[3],h[3],0,0,h[1]]; Plotter:-new_color(1,0.8,0);Plotter:-new_thickness();Plotter:-plot_element(a,b,c) else x:=[seq(nods[n[i]][1]+h[i][1],i=1..3)]; y:=[seq(nods[n[i]][2]+h[i][2],i=1..3)]; Plotter:-new_color(rgb[]);Plotter:-new_thickness();Plotter:-plot_element(x,y) end if; end do; NULL end proc: animate_displacements:=proc(scale,llc::list,urc::list) global nods,disp,elems; local view,frms,e,i,j,k,l,x,y,mxv,scl,ssc,n,xlgth,ylgth,lgth; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); mxv:=max(seq(disp[i][1]^2+disp[i][2]^2,i=1..nops(disp))); ssc:=1;if nargs>0 then ssc:=max(ssc,args[1]) end if; scl:=evalf(lgth/sqrt(mxv))/ssc; frms:=50;Plotter:-start_2D_animation();Plotter:-draw_title(`Deformed Mesh`); for i to frms do j:=2.0*i/frms;if j>1 then j:=2-j end if; Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for e in elems do n:=e[1..3]; x,y:=seq([seq(nods[n[k]][l],k=1..3)],l=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[k],y[k]],1),k=1..3) end do; [seq([seq(disp[k][l]*scl*j,l=1..2)],k=1..nops(nods))],[1,0.7,0]; deform_mesh(%); Plotter:-upgrade_frames() end do; if nargs>1 then view:=args[2..3] else view:=NULL end if; Plotter:-show_2D_animation(view); end proc: plot_displacements:=proc(scale,llc::list,urc::list) global nods,disp; local elem,x,y,mxv,mnv,dv,scl,ssc,n,v,i,j,lab,view,xlgth,ylgth,lgth,th,a; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); ssc:=1;if nargs>0 then ssc:=max(ssc,args[1]) end if; lgth:=max(xlgth,ylgth)/ssc; Plotter:-start_2D_plot();Plotter:-draw_title(`Displacement Field`); Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for elem in elems do n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; mxv:=max(seq(disp[i][1]^2+disp[i][2]^2,i=1..nops(disp))); mxv:=sqrt(mxv);scl:=evalf(lgth/mxv); mnv:=min(seq(disp[i][1]^2+disp[i][2]^2,i=1..nops(disp))); mnv:=sqrt(mnv)*scl;dv:=evalf(abs(lgth-mnv)); for i to nops(nods) do n:=nods[i];v:=[disp[i][1]*scl,disp[i][2]*scl]; th:=sqrt(v[1]^2+v[2]^2); if dv<10^(1-Digits) then a:=-1 else a:=evalf(cos(abs(th-mnv)/dv*Pi)) end if; Plotter:-new_color(rgb_color(a));Plotter:-new_arrow(n,v,th*0.25,th*0.5,1/3) end do; lab:=cat("Maximum Value: ",convert(evalf(mxv,2),string)); if nargs>1 then view:=args[2..3] else view:=NULL;Plotter:-draw_label([lab],[[1,0,0]]) end if; Plotter:-show_2D_plot(view); end proc: rgb_color:=proc(a) local x,n,m,r,g,b; x:=min(a,1);x:=max(x,-1); n:=1.3;m:=1.5; r:=1/(1+n**2*(x+1)**2); g:=1/(1+m**2*x**2); b:=1/(1+n**2*(x-1)**2); r,g,b end proc: plot_stresses:=proc(strs::list,scale,llc::list,urc::list) global nods,elems; local nval,n,x,y,xg,yg,mxv,scl,ssc,k,e,i,j,lab,view,xlgth,ylgth,lgth; if nops(strs)=nops(nods) then nval:=true else nval:=false end if; xlgth:=max(seq(nods[i][1],i=1..nops(nods))); xlgth:=0.125*(xlgth-min(seq(nods[i][1],i=1..nops(nods)))); ylgth:=max(seq(nods[i][2],i=1..nops(nods))); ylgth:=0.125*(ylgth-min(seq(nods[i][2],i=1..nops(nods)))); lgth:=max(xlgth,ylgth); mxv:=max(seq(seq(abs(strs[i][j]),j=1..2),i=1..nops(strs))); ssc:=1;if nargs>1 then ssc:=max(ssc,args[2]) end if; scl:=evalf(lgth/mxv)/ssc; Plotter:-start_2D_plot();Plotter:-draw_title(convert(readstat("plot title: "),string)); Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for e in elems do Plotter:-new_color();n:=e[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; k:=0; if nval then for n in nods do Plotter:-new_color();k:=k+1; stress_arrow(n,[strs[k][1]*scl,strs[k][2]*scl,strs[k][3]*Pi/180]) end do else for e in elems do Plotter:-new_color();k:=k+1;n:=e[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); xg,yg:=add(x[i],i=1..3)/3,add(y[i],i=1..3)/3; stress_arrow([xg,yg],[strs[k][1]*scl,strs[k][2]*scl,strs[k][3]*Pi/180]) end do end if; lab:=cat("Maximum Value: ",convert(evalf(mxv,2),string)); if nargs>2 then view:=args[3..4] else view:=NULL;Plotter:-draw_label([lab],[[0,0,1]]) end if; Plotter:-show_2D_plot(view); end proc: stress_arrow:=proc(n::list,si::list) global nods; local uu,j,u,v,th,t1,t2; uu:=[evalf(si[1]),evalf(si[2]),[evalf(si[3]),evalf(si[3]+Pi/2)]]; for j to 2 do u,v:=uu[j]*cos(uu[3][j]),uu[j]*sin(uu[3][j]); th:=sqrt(u^2+v^2); if j=1 then Plotter:-new_color(1,0,0);t1,t2:=th*0.8,1/3 else Plotter:-new_color(0,1,1);t1,t2:=0,0 end if; Plotter:-new_arrow(n,[u,v],th/5,t1,t2); Plotter:-new_arrow(n,[-u,-v],th/5,t1,t2) end do; NULL end proc: plot_contour_lines:=proc(nlevels::list,lvls::list,llc::list,urc::list) global elems,nods; local elem,level,levels,line,view,x,y,p_in,p_out,d,dm,n,z,i,j,k,cols,a,d0; Plotter:-start_2D_plot(); Plotter:-draw_title(convert(readstat("plot title: "),string)); if nargs>1 then levels:=lvls else d:=evalf((max(nlevels[])-min(nlevels[]))/6); levels:=[seq(min(nlevels[])+i*d,i=1..5)] end if; Plotter:-new_color();Plotter:-new_pen_style();Plotter:-new_thickness();Plotter:-new_axes(0); for elem in elems do; n:=elem[1..3]; x,y:=seq([seq(nods[n[i]][j],i=1..3)],j=1..2); Plotter:-move_pen([x[3],y[3]],0);seq(Plotter:-move_pen([x[i],y[i]],1),i=1..3) end do; Plotter:-new_pen_style();Plotter:-new_thickness(3); d0:=min(levels[]);dm:=max(levels[])-d0;cols:=[]; if abs(dm)<10000*10^(1-Digits) then RETURN(`identical levels`) end if; for level in levels do d:=(level-d0)/dm;a:=evalf(cos(d*Pi)); Plotter:-new_color(rgb_color(a));cols:=[cols[],[rgb_color(a)]]; line:=[]; for elem in elems do n:=elem[1..3];z:=seq(nlevels[n[i]],i=1..3); for i to 3 do j,k:=loop_123(i); if z[i]<=level and z[j]>=level then p_in:=crossing_point([n[j],n[i]],[z[j],z[i]],level); if z[k][] then for i to nops(line) do Plotter:-move_pen(line[i][1],0);Plotter:-move_pen(line[i][2],1) end do end if end do; if nargs>2 then view:=args[3..4] else view:=NULL; [seq(convert(evalf(levels[i],2),string),i=1..nops(levels))],cols; Plotter:-draw_label(%) end if; Plotter:-show_2D_plot(view); end proc: crossing_point:=proc(n::list,z::list,level::algebraic) global nods; local t,i; if abs(z[1]-z[2])<10^(1-Digits) then t:=[1] else t:=[evalf((level-z[2])/(z[1]-z[2]))] end if; t:=[t[],evalf(1-t[1])]; [add(t[i]*nods[n[i]][1],i=1..2),add(t[i]*nods[n[i]][2],i=1..2)] end proc: loop_123:=proc(i::posint) local j,k; if i<3 then j:=i+1 else j:=1 end if; if j<3 then k:=j+1 else k:=1 end if; j,k end proc: end module: savelibname := "."; savelib( 'Plotter' ); savelib( 'Cgt_fem' ); savelib( 'G_cgt_fem' ); savelib( 'Cst_fem' ); savelib( 'G_cst_fem' ); libname:=savelibname, libname;