function lambda = c0P2(h)
%% Mesh
[node,elem] = squarequadmesh([,,,],h);
elem = elem(:,[,,,]);
showmesh(node,elem);
findnode(node);
findquadelem(node,elem);
%% Construct Data Structure
[elem2dof,edge,inDof] = c0dofP2(elem);
elem2dof=double(elem2dof);
N = size(node,); NT = size(elem,);
Ndof = N+NT+size(edge,);
A=sparse(Ndof,Ndof);
B=sparse(Ndof,Ndof);
%% Assemble stiffness matrix
% Since Dphi_i*Dphi_j is quadratic,
% numerical quadrature rule is used here
option.quadorder = ; % default order
[pts,weight] = quadquadpts(option.quadorder);
pts=pts*-;
x=pts(:,);y=pts(:,);
h1=h/;h2=h1;
area=h2*h1;
%% Assemble Matrix
for i=:
for j=i:
DuDv=area*(Dxphi(x,y,h1,i).*Dxphi(x,y,h1,j)...
+Dyphi(x,y,h2,i).*Dyphi(x,y,h2,j))'*weight;
uv=area*(phi(x,y,i).*phi(x,y,j))'*weight;
if i==j
A = A + sparse(elem2dof(:,i),elem2dof(:,j),DuDv,Ndof,Ndof);
B = B + sparse(elem2dof(:,i),elem2dof(:,j),uv,Ndof,Ndof);
else
A = A + sparse([elem2dof(:,i);elem2dof(:,j)],...
[elem2dof(:,j);elem2dof(:,i)],DuDv,Ndof,Ndof);
B = B + sparse([elem2dof(:,i);elem2dof(:,j)],...
[elem2dof(:,j);elem2dof(:,i)],uv,Ndof,Ndof);
end
end
end
%% Solve Ax = lambda Bx, and its first solution is *pi^
lambda=eigs(A(inDof,inDof),B(inDof,inDof),,'sm');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dxphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = Dxphi(xi,eta,L1,i) % gradient of basis phi for x
switch i
case
s = (eta.*(.*xi - ).*(eta - ))/(*L1);
case
s = (eta.*(.*xi + ).*(eta - ))/(*L1);
case
s = (eta.*(.*xi + ).*(eta + ))/(*L1);
case
s = (eta.*(.*xi - ).*(eta + ))/(*L1);
case
s = -(eta.*xi.*(eta - ))/L1;
case
s = -((eta.^ - ).*(*xi + ))/(*L1);
case
s = -(eta.*xi.*(eta + ))/L1;
case
s = -((eta.^ - ).*(.*xi - ))/(*L1);
case
s = (*xi.*(eta.^ - ))/L1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dxphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = Dyphi(xi,eta,L2,i) % gradient of basis phi for y
switch i
case
s = (xi.*(.*eta - ).*(xi - ))/(*L2);
case
s = (xi.*(.*eta - ).*(xi + ))/(*L2);
case
s = (xi.*(.*eta + ).*(xi + ))/(*L2);
case
s = (xi.*(.*eta + ).*(xi - ))/(*L2);
case
s = -((.*eta - ).*(xi.^ - ))/(*L2);
case
s = -(eta.*xi.*(xi + ))/L2;
case
s =-((*eta + ).*(xi.^ - ))/(*L2);
case
s = -(eta.*xi.*(xi - ))/L2;
case
s = (*eta.*(xi.^ - ))/L2;
end
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction phi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = phi(xi,eta,i) % gradient of basis phi
switch i
case
s = (eta.*xi.*(eta - ).*(xi - ))/;
case
s = (eta.*xi.*(eta - ).*(xi + ))/;
case
s = (eta.*xi.*(eta + ).*(xi + ))/;
case
s = (eta.*xi.*(eta + ).*(xi - ))/;
case
s = -(eta.*(xi.^ - ).*(eta - ))/;
case
s = -(xi.*(eta.^ - ).*(xi + ))/;
case
s = -(eta.*(xi.^ - ).*(eta + ))/;
case
s = -(xi.*(eta.^ - ).*(xi - ))/;
case
s = (eta.^ - ).*(xi.^ - );
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
lambda = c0P2(h)
function [elem2dof,edge,inDof] = c0dofP2(elem)
totalEdge=sort([elem(:,[,]);elem(:,[,]);elem(:,[,]);elem(:,[,])],);
[edge, i2, j] = myunique(totalEdge);
N = max(elem(:));
NT = size(elem,);
NE = size(edge,);
elem2edge = reshape(j,NT,);
elem2dof = uint32([elem N+elem2edge (N+NE+:N+NE+NT)']);
i1(j(*NT:-:)) = *NT:-:;
i1 = i1';
bdEdgeIdx = (i1 == i2);
isBdDof = false(N+NE+NT,);
isBdDof(edge(bdEdgeIdx,:)) = true; % nodal
idx = find(bdEdgeIdx);
isBdDof(N+idx) = true;
inDof = find(~isBdDof);
end
[elem2dof,edge,inDof] = c0dofP2(elem)
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