EFG.zip_EFG_cantilever beam_run

%Created by J.T.B. Overvelde - 18 April 2012 %Master's thesis - The Moving Noda Approach in Topology Optimization %http://www.overvelde.com % %Initialize problem, nodal distribution and background mesh. function time1=InitMesh() GlobalConst tic %Mesh timer %Problem definition constants pCon.Lx=48; %Width of the beam pCon.Ly=12; %Height of the beam pCon.P=0.1; %Maximum force in parabolic loading pCon.E=1; %Elasticiy modulus pCon.nu=0.3; %Poisson's ratio pCon.b=[0;0]; %Body Force pCon.D=pCon.E/(1-pCon.nu^2)*[1 pCon.nu 0; %Plain stress elasticy matrix pCon.nu 1 0; 0 0 (1-pCon.nu)/2]; %Meshless discretization constants mCon.nx=10; %number of nodes along the width mCon.ny=10; %number of nodes along the height mCon.d=2.5; %Relative smoothing length mCon.pn=2; %Number of monomials in the Monomial basis vector mCon.nG=2; %Number of quadrature points in one dimension in each cells mCon.BCuType=2; %Displacement Boundary type, 1: integration cells 2: Collocated method mCon.DispNodesType=1; %Nodal distribution type, 1: pattern 2: random mCon.n=mCon.nx*mCon.ny; %Number of total nodes mCon.dx=pCon.Lx/(mCon.nx-1); %Horizontal distance between nodes mCon.dy=pCon.Ly/(mCon.ny-1); %Vertical distance between mCon.mx=mCon.nx-1; %Number of cells along the width mCon.my=mCon.ny-1; %Number of cells along the height mCon.m=mCon.mx*mCon.my; %number of total internal cells if mCon.BCuType==1 mCon.mbl=mCon.my; %Number of boundary cells along the left edge mCon.mblc=0; %Number of collocated boundary particles along left edge mCon.me=mCon.mbl+1; %Number of essential boundary condition nodes else mCon.mbl=0; %Number of boundary cells along the left edge mCon.mblc=mCon.ny; %Number of collocated boundary particles along left edge mCon.me=0; %Number of essential boundary condition nodes end mCon.mbr=mCon.my; %Number of boundary cells along the right edge mCon.mb=mCon.mbl+mCon.mbr; %Number of total boundary cells. mCon.dm=[mCon.d*mCon.dx ; mCon.d*mCon.dy]; %Smoothing length in x and y direction, respectively. mCon.drn=1; %Position in domain for which the final solution is calculated sCon.snx=25; %Number of nodes in x direction for which the solution is calculated sCon.sny=25; %Number of nodes in y direction for which the solution is calculated sCon.sn=sCon.snx*sCon.sny; %Number of nodes for which the solution is calculated sCon.dsx=pCon.Lx/(sCon.snx-1); %Horizontal distance between solution nodes sCon.dsy=pCon.Ly/(sCon.sny-1); %Vertical distance between solution nodes %Create nodes for i=1:mCon.n nodes(i).x=[]; %nodes coordinates nodes(i).nen=[]; %Neighboring nodes end %Create Solution nodes for i=1:sCon.sn solnodes(i).x=[]; %solnodes coordinates solnodes(i).nen=[]; %Neighboring nodes solnodes(i).nx=[]; %Matrix x location solnodes(i).ny=[]; %Matrix y location end %Create internal background cells for i=1:mCon.m cells(i).x=[]; %coordinates of cells center cells(i).dx=[]; %width and height cells center cells(i).ni=mCon.nG*mCon.nG; %number of integration points cells(i).J=[]; %Jacobian for j=1:cells(i).ni cells(i).int(j).x=[]; %Coordinates of integration points cells(i).int(j).nen=[]; %Neighboring nodes cells(i).int(j).w=[]; %Weight values for each intergration point cells(i).int(j).cv=[]; %Body force vector end end %Create left and right boundary cells for i=1:mCon.mb bcells(i).BC=0; %Boundary type - 0 is undefined bcells(i).x=[]; %cells center bcells(i).dx=[]; %width and height cells center bcells(i).ni=mCon.nG; %number of integration points bcells(i).J=[]; %Jacobian bcells(i).Ni=[]; %Two neighboring element numbers - only necessary for essential boundary condition for j=1:bcells(i).ni bcells(i).int(j).x=[]; %Integration points bcells(i).int(j).nen=[]; %Neighboring nodes bcells(i).int(j).w=[]; %Weight values for each intergration point bcells(i).int(j).bv=[]; %Boundary specific scalar or vector bcells(i).int(j).N=zeros(2,1); %Interpolation value - only necessary for essential boundary condition end end ag=1; %Constant telling if the material distribution is homogeneous while ag==1 ag=0; if mCon.DispNodesType==1 %Organized distribution of particles for i=1:mCon.nx for j=1:mCon.ny nodes((i-1)*mCon.ny+j).x=[mCon.dx*(i-1); mCon.dy*(j-1)-pCon.Ly/2]; end end elseif mCon.DispNodesType==2 %Random distribution in domain for i=1:mCon.n nodes(i).x=[pCon.Lx*rand(1); pCon.Ly*rand(1)-pCon.Ly/2]; end elseif mCon.DispNodesType==3 %First Organize the nodes, then distribute them randomly with small variation around node for i=1:mCon.nx for j=1:mCon.ny nodes((i-1)*mCon.ny+j).x=[mCon.dx*(i-1); mCon.dy*(j-1)-pCon.Ly/2]; end end for i=1:mCon.n a=1; while a==1 a=0; dx=mCon.drn*[mCon.dx*(rand(1)-0.5); mCon.dy*(rand(1)-0.5)]; if nodes(i).x(2)+dx(2)>pCon.Ly/2 a=1; elseif nodes(i).x(2)+dx(2)<-pCon.Ly/2 a=1; else nodes(i).x=nodes(i).x+dx; end end end end %solnodes coordinates for i=1:sCon.snx for j=1:sCon.sny solnodes((i-1)*sCon.sny+j).x=[sCon.dsx*(i-1); sCon.dsy*(j-1)-pCon.Ly/2]; solnodes((i-1)*sCon.sny+j).nx=i; solnodes((i-1)*sCon.sny+j).ny=j; end end %Boundary type for i=1:mCon.mb if i<=mCon.mbl bcells(i).BC=1; %Essential boundary condition else bcells(i).BC=2; %Traction boundary condition end end [t,w]=ConGauss(mCon.nG); %Gaus quadrature in general coordinates %Internal cells coordinates, integration points and weights for i=1:mCon.mx for j=1:mCon.my cells((i-1)*mCon.my+j).x=[mCon.dx/2+mCon.dx*(i-1); mCon.dy/2+mCon.dy*(j-1)-pCon.Ly/2]; cells((i-1)*mCon.my+j).dx=[mCon.dx;mCon.dy]; x1=cells((i-1)*mCon.my+j).x(1)-cells((i-1)*mCon.my+j).dx(1)/2; x2=cells((i-1)*mCon.my+j).x(1)+cells((i-1)*mCon.my+j).dx(1)/2; y1=cells((i-1)*mCon.my+j).x(2)-cells((i-1)*mCon.my+j).dx(2)/2; y2=cells((i-1)*mCon.my+j).x(2)+cells((i-1)*mCon.my+j).dx(2)/2; Gcx=1/2*(x2+x1); Gmx=1/2*(x2-x1); Gcy=1/2*(y2+y1); Gmy=1/2*(y2-y1); cells((i-1)*mCon.my+j).J=Gmx*Gmy; for k=1:mCon.nG for l=1:mCon.nG cells((i