kfpro.m

用无网格法计算VAN DE POL 振子

function [K_r,K_o,F_r,F_o,PHI] = KFPRO(m,nnodes,xi,npoints,x,dm,wtype,para,disppen,curpen,delta_t,time,...
                                        yi, yh3, yh2, yh1, qi, qh2, qh1)
% KFRO IS A FUNCTION TO CALCULATE THE TOTAL STIFFNESS AND LOADS OF RISER AND OSSILATOR 
% BY USING MLS APPROXIMATION
%
% INPUT PARAMETERS
%    m - Total number of basis functions (1: Constant basis;  2: Linear basis;  3: Quadratic basis)
%    nnodes  - Total number of nodes used to construct MLS approximation
%    npoints - Total number of points whose MLS shape function to be evaluated
%    xi(nnodes) - Coordinates of nodes used to construct MLS approximation
%    x(npoints) - Coordinates of points whose MLS shape function to be evaluated
%    dm(nnodes) - Radius of support of nodes
%    wtype - Type of weight function
%    para  - Weight function parameter
%
% OUTPUT PARAMETERS
%    PHI   - MLS Shpae function
%    DPHI  - First order derivatives of MLS Shpae function
%    D2PHI - Second order derivatives of MLS Shpae function
%
% INITIALIZE WEIGHT FUNCTION MATRICES
wi   = zeros (1, nnodes);  % Weight funciton
dwi  = zeros (1, nnodes);
d2wi = zeros (1, nnodes);

% INITIALIZE SHAPE FUNCTION MATRICES
PHI   = zeros(npoints, nnodes);
DPHI  = zeros(npoints, nnodes);
D2PHI = zeros(npoints, nnodes);

K_r   = zeros(npoints, nnodes);
K_o   = zeros(npoints, nnodes);
F_r   = zeros(npoints, 1);
F_o   = zeros(npoints, 1);

% LOOP OVER ALL EVALUATION POINTS TO CALCULATE VALUE OF SHAPE FUNCTION Fi(X)
for j = 1 : npoints

	% DETERMINE WEIGHT FUNCTIONS AND THEIR DERIVATIVES AT EVERY NODE
	for i = 1 : nnodes
		di = x(j) - xi(i);
      [wi(i),dwi(i),d2wi(i),d3wi(i)] = Weight(wtype, para, di, dm(i));
	end
   
   % EVALUATE BASIS p, B MATRIX AND THEIR DERIVATIVES
   
   if (m == 3)
      p = [ones(1, nnodes); xi; xi.^2]; 
      px   = [1; x(j); x(j).^2];
      dpx  = [0; 1; 2*x(j)];
      d2px = [0; 0; 2];
      d3px = [0; 0; 0];
%//       d4px = [0; 0; 0];
      
      B    = p .* [wi; wi; wi];
      DB   = p .* [dwi; dwi; dwi];
      D2B  = p .* [d2wi; d2wi; d2wi];
      D3B  = p .* [d3wi; d3wi; d3wi];
%//       D4B  = p .* [d4wi; d4wi; d4wi];
   else
      error('Invalid order of basis.');
   end
   
   % EVALUATE MATRICES A AND ITS DERIVATIVES
	A   = zeros (m, m);
	DA  = zeros (m, m);
	D2A = zeros (m, m);
    D3A = zeros (m, m);
%//     D4A = zeros (m, m);
	for i = 1 : nnodes
      pp = p(:,i) * p(:,i).';
      
      A   = A   + wi(i) * pp;
      DA  = DA  + dwi(i) * pp;
      D2A = D2A + d2wi(i) * pp;
      D3A = D3A + d3wi(i) * pp;
%//       D4A = D4A + d4wi(i) * pp;
    end
   
   AInv = inv(A);
      
   rx  = AInv * px;
   PHI(j,:) = rx.' * B;   % shape function
   drx  = AInv * (dpx -DA * rx);
   DPHI(j,:) = drx.' * B + rx.' * DB;                % first order derivatives of shape function
   
   d2rx  = AInv * (d2px - 2 * DA * drx - D2A * rx);
   D2PHI(j,:) = d2rx.'*B + 2*drx.'*DB + rx.'*D2B;             % second order derivatives of shape function
   
   d3rx  = AInv * (d3px - 3*D2A*drx - 3*DA*d2rx - D3A*rx);
   D3PHI(j,:) = d3rx.'*B + 3*d2rx.'*DB + 3*drx.'*D2B + rx.' * D3B;     % third order derivatives of shape function ????????????
   
%//    d4rx  = AInv * (d4px - 4*D3A*drx - 6*D2A*d2rx -4*DA*d3rx - D4A*rx);
%//    D4PHI(j,:) = d4rx.'*B + 4*d3rx.'*DB + 6*d2rx.'*D2B + 4*drx.'*D3B + rx.'*D4B; % fourth order derivatives of shape function
   
   % ACCESS THE TOTAL DAMPING, MASS AND SOME OTHER CONSTANTS
   [m_i,v_i,EI,eps,ome,Amp,D_out,const1,const2,C_t,D_t,m_t] = nodes_CC(yi,time,x,npoints); 
   
   K1   = m_t/delta_t^2.*PHI(j,:) + C_t/delta_t.*PHI(j,:) + 2*m_i*v_i/delta_t.*DPHI(j,:);
   D2K1 = m_t/delta_t^2.*D2PHI(j,:) + C_t/delta_t.*D2PHI(j,:) + 2*m_i*v_i/delta_t.*D3PHI(j,:);
   K2   = 1/delta_t^2.*PHI(j,:) + eps/delta_t.*ome.*(qi.'.^2-ones(1,npoints)).*PHI(j,:);
   % ASSEMBLED TO STIFFNESS MATRICES
   K_r = K_r + K1.'*K1;
   K_o = K_o + K2.'*K2;
   
   F1yh2  = 2*m_t/delta_t^2.*PHI(j,:) + C_t/delta_t.*PHI(j,:) + (2*m_i*v_i/delta_t+const2).*DPHI(j,:) - const1.*D2PHI(j,:);
   F1yh2_ = EI .* D2PHI(j,:);
   F1yh1  = - m_t/delta_t^2.*PHI(j,:);
   F1qh2  = D_t.*PHI(j,:);
   
   F2qh2  = eps/delta_t.*ome.*(qi.'.^2).*PHI(j,:) + (2/delta_t^2.*ones(1,nnodes)-eps/delta_t.*ome-ome.^2).*PHI(j,:);
   F2qh1  = - 1/delta_t^2.*PHI(j,:);
   F2D2y  = Amp/D_out/delta_t^2.*PHI(j,:);
   % ASSEMBLED TO LOAD VECTORS
   F_r = F_r + (K1.' * F1yh2 + D2K1.' * F1yh2_)*yh2 + (K1.' * F1yh1)*yh1 + (K1.' * F1qh2)*qh2;
   F_o = F_o + (K2.' * F2qh2)*qh2 + (K2.' * F2qh1)*qh1 + (K2.' * F2D2y)*(yh3-2*yh2-yh1);

end

% ENFORCEMENT OF DISPLACEMENT BOUDNARY CONDITIONS
K_r = K_r + disppen*(PHI(1,:).'*PHI(1,:) + PHI(npoints,:).'*PHI(npoints,:)) + ...
      curpen*(D2PHI(1,:).'*D2PHI(1,:) + D2PHI(npoints,:).'*D2PHI(npoints,:));
K_o = K_o + disppen*(PHI(1,:).'*PHI(1,:) + PHI(npoints,:).'*PHI(npoints,:)) + ...
      curpen*(D2PHI(1,:).'*D2PHI(1,:) + D2PHI(npoints,:).'*D2PHI(npoints,:));