EFG_beam.rar_EFG_efg matlab code_matlab efg

% +-----------------------------------------------------------------------+ % |-------------- EFG Meshless for 1D Bending beam problem --------------| % |------- Copyright: Canh V. Le - Harm Askes - Matthew Gilbert ----------| % |------------------- The University of Sheffield 2007 ----------------- | % +-----------------------------------------------------------------------+ clear all % +------------+ % | Input data | % +------------+ length = 10; % The length of bar (m) E = 1000; % Young Modulus I = 1000; q = 20; % Load % Parameter for exact solution u0 = 0; theta0 = 0; M0 = 0; Q0 = 0; % +------------------+ % | Number of nodes: | % +------------------+ nno = 41; for i = 1:nno xno(i) = (i-1)*length/(nno-1); % Coordinate of nodes vno(i) = length/(nno-1); % Distance between two nodes end % +------------------------------------------------+ % | Number of integration points (Gaussian point): | % +------------------------------------------------+ nip = nno*2; % The number of Gaussian points for i=1:nip xip(i) = (i-0.5)*length/nip; % Coordinate of Gaussian points vip(i) = length/nip; % Distance between two Gaussian ponits end % +--------------------------------+ % | Radius of domain of influence: | % +--------------------------------+ dinflu = 3*(xno(2)-xno(1)); % vno=x2-x1 % +----------------------------------------------------+ % | Constructing shape functions (at Gaussian points): | % +----------------------------------------------------+ [phi0,phi1,phi2] = shapef_beam_LU(xno,xip,dinflu); % storing | phi_1(x_1) phi_1(x_2) ... phi_1(x_Gauss) | % | phi_2(x_1) phi_2(x_2) ... phi_2(x_Gauss) | % ... % | phi_n(x_1) phi_n(x_2) ... phi_n(x_Gauss) | n =nno % +--------------------------------------+ % | Essential boundary condition at x=0: | % +--------------------------------------+ [phi0e,phi1e,phi2e] = shapef_beam_LU(xno,0,dinflu); % +------------------------------------+ % | Natural boundary condition at x=L: | % +------------------------------------+ [phi0n,phi1n,phi2n] = shapef_beam_LU(xno,length,dinflu); % +---------------------------+ % | Global stiffnes matrix K: | % +---------------------------+ K = stiffn_beam(E,I,phi2,phi0e,phi1e,vip,nno,nip); % +------------------------+ % | Global force vector f: | % +------------------------+ f = rhs_beam(phi0,phi0n,phi1n,q,Q0,M0,u0,theta0,vip,nno,nip); % +-------------------+ % | Solving equation: | % +-------------------+ d = inv(K) * f; % +------------------------+ % | Output representation: | % +------------------------+ % Plot result: displacement at Gaussian points figure; h0 = gcf; set(h0,'name','EFG Displacement'); set(h0,'NumberTitle','off'); dipw = phi0' * (d(1:nno)); plot (xip,dipw,'o r') hold on %Plot exact solution: [uy,theta] = uexact_beam(xip,q,E,I,M0,Q0,length,u0,theta0); plot (xip,uy,'b'); figure; h1 = gcf; set(h1,'name','EFG Rotation'); set(h1,'NumberTitle','off'); dipd = phi1' * (d(1:nno)); plot (xip,dipd,'o r') hold on plot (xip,theta,'b') L2w=sqrt((uy-dipw')*(uy-dipw')') L2d=sqrt((uy-dipd')*(uy-dipd')') % +-----------------------------------------------------------------------+ % | --------------------------------The end------------------------------ | % +-----------------------------------------------------------------------+