https://github.com/team-pancho/HDG3D
Tip revision: af0d6cb24c63ad6fa07331f26a5733116db800b7 authored by shukaidu on 29 April 2020, 21:46:21 UTC
Merge pull request #2 from team-pancho/hdg3d_skdu
Merge pull request #2 from team-pancho/hdg3d_skdu
Tip revision: af0d6cb
Stiffness.m
function S=Stiffness(T,k,formula)
% S=stiffness(T,k,formula)
%
% Input:
% T : tetrahedrization
% k : polynomial degree
% formula : quadrature formula
% Output:
% S : d3 x d3 x Nelts stiffness matrix
% Last modified: March 27, 2013
% Computations in the reference element
xhat=formula(:,2);
yhat=formula(:,3);
zhat=formula(:,4);
[~,px,py,pz]=dubiner3d(2*xhat-1,2*yhat-1,2*zhat-1,k);
px=2*px; wpx=bsxfun(@times,formula(:,5),px);
py=2*py; wpy=bsxfun(@times,formula(:,5),py);
pz=2*pz; wpz=bsxfun(@times,formula(:,5),pz);
Kxx=(1/6)*px'*wpx;
Kxy=(1/6)*px'*wpy;
Kxz=(1/6)*px'*wpz;
Kyy=(1/6)*py'*wpy;
Kyz=(1/6)*py'*wpz;
Kzz=(1/6)*pz'*wpz;
% Geometric constructions
x12=T.coordinates(T.elements(:,2),1)-T.coordinates(T.elements(:,1),1);
y12=T.coordinates(T.elements(:,2),2)-T.coordinates(T.elements(:,1),2);
z12=T.coordinates(T.elements(:,2),3)-T.coordinates(T.elements(:,1),3);
x13=T.coordinates(T.elements(:,3),1)-T.coordinates(T.elements(:,1),1);
y13=T.coordinates(T.elements(:,3),2)-T.coordinates(T.elements(:,1),2);
z13=T.coordinates(T.elements(:,3),3)-T.coordinates(T.elements(:,1),3);
x14=T.coordinates(T.elements(:,4),1)-T.coordinates(T.elements(:,1),1);
y14=T.coordinates(T.elements(:,4),2)-T.coordinates(T.elements(:,1),2);
z14=T.coordinates(T.elements(:,4),3)-T.coordinates(T.elements(:,1),3);
x12=x12'; x13=x13'; x14=x14';
y12=y12'; y13=y13'; y14=y14';
z12=z12'; z13=z13'; z14=z14'; % row vectors with differences of coords
idetB=1./(x12.*y13.*z14-x12.*y14.*z13-y12.*x13.*z14...
+y12.*x14.*z13+z12.*x13.*y14-z12.*x14.*y13);
dxdx=kron(idetB.*(y13.*z14-z13.*y14).^2,Kxx)...
+kron(idetB.*(y13.*z14-z13.*y14).*(-y12.*z14+z12.*y14),Kxy+Kxy')...
+kron(idetB.*(y13.*z14-z13.*y14).*(y12.*z13-z12.*y13),Kxz+Kxz')...
+kron(idetB.*(-y12.*z14+z12.*y14).^2,Kyy)...
+kron(idetB.*(-y12.*z14+z12.*y14).*(y12.*z13-z12.*y13),Kyz+Kyz')...
+kron(idetB.*(y12.*z13-z12.*y13).^2,Kzz);
dydy=kron(idetB.*(-x13.*z14+z13.*x14).^2,Kxx)...
+kron(idetB.*(-x13.*z14+z13.*x14).*(x12.*z14-z12.*x14),Kxy+Kxy')...
+kron(idetB.*(-x13.*z14+z13.*x14).*(-x12.*z13+z12.*x13),Kxz+Kxz')...
+kron(idetB.*(x12.*z14-z12.*x14).^2,Kyy)...
+kron(idetB.*(x12.*z14-z12.*x14).*(-x12.*z13+z12.*x13),Kyz+Kyz')...
+kron(idetB.*(-x12.*z13+z12.*x13).^2,Kzz);
dzdz=kron(idetB.*(x13.*y14-y13.*x14).^2,Kxx)...
+kron(idetB.*(x13.*y14-y13.*x14).*(-x12.*y14+y12.*x14),Kxy+Kxy')...
+kron(idetB.*(x13.*y14-y13.*x14).*(x12.*y13-y12.*x13),Kxz+Kxz')...
+kron(idetB.*(-x12.*y14+y12.*x14).^2,Kyy)...
+kron(idetB.*(-x12.*y14+y12.*x14).*(x12.*y13-y12.*x13),Kyz+Kyz')...
+kron(idetB.*(x12.*y13-y12.*x13).^2,Kzz);
S=dxdx+dydy+dzdz;
d3=nchoosek(k+3,3);
Nelts=size(T.elements,1);
S=reshape(S,[d3 d3 Nelts]);
return
% These matrices are not needed
dxdy=kron(idetB.*(y13.*z14-z13.*y14).*(-x13.*z14+z13.*x14),Kxx)...
+kron(idetB.*(y13.*z14-z13.*y14).*(x12.*z14-z12.*x14),Kxy)...
+kron(idetB.*(y13.*z14-z13.*y14).*(-x12.*z13+z12.*x13),Kxz)...
+kron(idetB.*(-x13.*z14+z13.*x14).*(-y12.*z14+z12.*y14),Kxy')...
+kron(idetB.*(-y12.*z14+z12.*y14).*(x12.*z14-z12.*x14),Kyy)...
+kron(idetB.*(-y12.*z14+z12.*y14).*(-x12.*z13+z12.*x13),Kyz)...
+kron(idetB.*(-x13.*z14+z13.*x14).*(y12.*z13-z12.*y13),Kxz')...
+kron(idetB.*(x12.*z14-z12.*x14).*(y12.*z13-z12.*y13),Kyz')...
+kron(idetB.*(y12.*z13-z12.*y13).*(-x12.*z13+z12.*x13),Kzz);
dxdz=kron(idetB.*(y13.*z14-z13.*y14).*(x13.*y14-y13.*x14),Kxx)...
+kron(idetB.*(y13.*z14-z13.*y14).*(-x12.*y14+y12.*x14),Kxy)...
+kron(idetB.*(y13.*z14-z13.*y14).*(x12.*y13-y12.*x13),Kxz)...
+kron(idetB.*(x13.*y14-y13.*x14).*(-y12.*z14+z12.*y14),Kxy')...
+kron(idetB.*(-y12.*z14+z12.*y14).*(-x12.*y14+y12.*x14),Kyy)...
+kron(idetB.*(-y12.*z14+z12.*y14).*(x12.*y13-y12.*x13),Kyz)...
+kron(idetB.*(x13.*y14-y13.*x14).*(y12.*z13-z12.*y13),Kxz')...
+kron(idetB.*(-x12.*y14+y12.*x14).*(y12.*z13-z12.*y13),Kyz')...
+kron(idetB.*(y12.*z13-z12.*y13).*(x12.*y13-y12.*x13),Kzz);
dydz=kron(idetB.*(-x13.*z14+z13.*x14).*(x13.*y14-y13.*x14),Kxx)...
+kron(idetB.*(-x13.*z14+z13.*x14).*(-x12.*y14+y12.*x14),Kxy)...
+kron(idetB.*(-x13.*z14+z13.*x14).*(x12.*y13-y12.*x13),Kxz)...
+kron(idetB.*(x13.*y14-y13.*x14).*(x12.*z14-z12.*x14),Kxy')...
+kron(idetB.*(x12.*z14-z12.*x14).*(-x12.*y14+y12.*x14),Kyy)...
+kron(idetB.*(x12.*z14-z12.*x14).*(x12.*y13-y12.*x13),Kyz)...
+kron(idetB.*(x13.*y14-y13.*x14).*(-x12.*z13+z12.*x13),Kxz')...
+kron(idetB.*(-x12.*y14+y12.*x14).*(-x12.*z13+z12.*x13),Kyz')...
+kron(idetB.*(-x12.*z13+z12.*x13).*(x12.*y13-y12.*x13),Kzz);