1 function lambda = c0P2(h) 2 %% Mesh 3 [node,elem] = squarequadmesh([0,1,0,1],h); 4 elem = elem(:,[1,4,3,2]); 5 showmesh(node,elem); 6 findnode(node); 7 findquadelem(node,elem); 8 %% Construct Data Structure 9 [elem2dof,edge,inDof] = c0dofP2(elem); 10 elem2dof=double(elem2dof); 11 N = size(node,1); NT = size(elem,1); 12 Ndof = N+NT+size(edge,1); 13 A=sparse(Ndof,Ndof); 14 B=sparse(Ndof,Ndof); 15 %% Assemble stiffness matrix 16 % Since Dphi_i*Dphi_j is quadratic, 17 % numerical quadrature rule is used here 18 option.quadorder = 4; % default order 19 [pts,weight] = quadquadpts(option.quadorder); 20 pts=pts*2-1; 21 x=pts(:,1);y=pts(:,2); 22 h1=h/2;h2=h1; 23 area=h2*h1; 24 %% Assemble Matrix 25 for i=1:9 26 for j=i:9 27 DuDv=area*(Dxphi(x,y,h1,i).*Dxphi(x,y,h1,j)... 28 +Dyphi(x,y,h2,i).*Dyphi(x,y,h2,j))'*weight; 29 uv=area*(phi(x,y,i).*phi(x,y,j))'*weight; 30 if i==j 31 A = A + sparse(elem2dof(:,i),elem2dof(:,j),DuDv,Ndof,Ndof); 32 B = B + sparse(elem2dof(:,i),elem2dof(:,j),uv,Ndof,Ndof); 33 else 34 A = A + sparse([elem2dof(:,i);elem2dof(:,j)],... 35 [elem2dof(:,j);elem2dof(:,i)],DuDv,Ndof,Ndof); 36 B = B + sparse([elem2dof(:,i);elem2dof(:,j)],... 37 [elem2dof(:,j);elem2dof(:,i)],uv,Ndof,Ndof); 38 end 39 end 40 end 41 %% Solve Ax = lambda Bx, and its first solution is 2*pi^2 42 lambda=eigs(A(inDof,inDof),B(inDof,inDof),1,'sm'); 43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 44 % subfunction Dxphi 45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 46 function s = Dxphi(xi,eta,L1,i) % gradient of basis phi for x 47 switch i 48 case 1 49 s = (eta.*(2.*xi - 1).*(eta - 1))/(4*L1); 50 case 2 51 s = (eta.*(2.*xi + 1).*(eta - 1))/(4*L1); 52 case 3 53 s = (eta.*(2.*xi + 1).*(eta + 1))/(4*L1); 54 case 4 55 s = (eta.*(2.*xi - 1).*(eta + 1))/(4*L1); 56 case 5 57 s = -(eta.*xi.*(eta - 1))/L1; 58 case 6 59 s = -((eta.^2 - 1).*(2*xi + 1))/(2*L1); 60 case 7 61 s = -(eta.*xi.*(eta + 1))/L1; 62 case 8 63 s = -((eta.^2 - 1).*(2.*xi - 1))/(2*L1); 64 case 9 65 s = (2*xi.*(eta.^2 - 1))/L1; 66 end 67 end 68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 69 % subfunction Dxphi 70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 71 function s = Dyphi(xi,eta,L2,i) % gradient of basis phi for y 72 switch i 73 case 1 74 s = (xi.*(2.*eta - 1).*(xi - 1))/(4*L2); 75 case 2 76 s = (xi.*(2.*eta - 1).*(xi + 1))/(4*L2); 77 case 3 78 s = (xi.*(2.*eta + 1).*(xi + 1))/(4*L2); 79 case 4 80 s = (xi.*(2.*eta + 1).*(xi - 1))/(4*L2); 81 case 5 82 s = -((2.*eta - 1).*(xi.^2 - 1))/(2*L2); 83 case 6 84 s = -(eta.*xi.*(xi + 1))/L2; 85 case 7 86 s =-((2*eta + 1).*(xi.^2 - 1))/(2*L2); 87 case 8 88 s = -(eta.*xi.*(xi - 1))/L2; 89 case 9 90 s = (2*eta.*(xi.^2 - 1))/L2; 91 end 92 end 93 94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 95 % subfunction phi 96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 97 function s = phi(xi,eta,i) % gradient of basis phi 98 switch i 99 case 1 100 s = (eta.*xi.*(eta - 1).*(xi - 1))/4; 101 case 2 102 s = (eta.*xi.*(eta - 1).*(xi + 1))/4; 103 case 3 104 s = (eta.*xi.*(eta + 1).*(xi + 1))/4; 105 case 4 106 s = (eta.*xi.*(eta + 1).*(xi - 1))/4; 107 case 5 108 s = -(eta.*(xi.^2 - 1).*(eta - 1))/2; 109 case 6 110 s = -(xi.*(eta.^2 - 1).*(xi + 1))/2; 111 case 7 112 s = -(eta.*(xi.^2 - 1).*(eta + 1))/2; 113 case 8 114 s = -(xi.*(eta.^2 - 1).*(xi - 1))/2; 115 case 9 116 s = (eta.^2 - 1).*(xi.^2 - 1); 117 end 118 end 119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 120 end
1 function [elem2dof,edge,inDof] = c0dofP2(elem) 2 totalEdge=sort([elem(:,[1,2]);elem(:,[2,3]);elem(:,[3,4]);elem(:,[4,1])],2); 3 [edge, i2, j] = myunique(totalEdge); 4 N = max(elem(:)); 5 NT = size(elem,1); 6 NE = size(edge,1); 7 elem2edge = reshape(j,NT,4); 8 elem2dof = uint32([elem N+elem2edge (N+NE+1:N+NE+NT)']); 9 i1(j(4*NT:-1:1)) = 4*NT:-1:1; 10 i1 = i1'; 11 bdEdgeIdx = (i1 == i2); 12 isBdDof = false(N+NE+NT,1); 13 isBdDof(edge(bdEdgeIdx,:)) = true; % nodal 14 idx = find(bdEdgeIdx); 15 isBdDof(N+idx) = true; 16 inDof = find(~isBdDof); 17 end
http://files.cnblogs.com/files/wangshixi12/%E5%8F%8C%E4%BA%8C%E6%AC%A1Lagrange%E6%9C%89%E9%99%90%E5%85%83.rar