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
projectHDG3d.m
function [Pqx,Pqy,Pqz,Pu]=projectHDG3d(T,coeffs,k,tau,formulas)
%[Pqx,Pqy,Pqz,Pu]=projectHDG3d(T,coeffs,k,tau,formulas)
%
%Input:
% T: Tetrahedrization
% coeffs: each cell is a vectorized function of three variables
% k: degree of polynomial
% formulas: quadrature matrix for mass matrix and convection matrix
%
%Output:
% Pqx: \Pi qx,
% Pqy: \Pi qy,
% Pqz: \Pi qz,
% Pu: \Pi u,
%Last modified: April 9, 2013
d2=nchoosek(k+2,2);
d3=nchoosek(k+3,3); %Nbasis
if k>0
d4=nchoosek(k+2,3); %dim P_{k-1}(K)
else
d4=0;
end
block3=@(x) (1+(x-1)*d3):(x*d3);
a=@(x,y,z) 1+0.*x;
Mass=MassMatrix(T,{a},k,formulas{1});
Mass=Mass{1};
rows=1+d4:d3;
Mass(rows,:,:)=[];
[tauPP,tauDP,nxDP,nyDP,nzDP]=matricesFace(T,tau,k,formulas{3});
Nelts=size(T.elements,1);
O=zeros(d4,d3,Nelts);
M=[Mass O O O ;
O Mass O O ;
O O Mass O ;
O O O Mass;
nxDP nyDP nzDP tauDP];
if k>0
Ints=testElem(coeffs,T,k-1,formulas{2});
qxP=Ints{1};
qyP=Ints{2};
qzP=Ints{3};
uP=Ints{4};
else
qxP=zeros(0,Nelts);
qyP=zeros(0,Nelts);
qzP=zeros(0,Nelts);
uP=zeros(0,Nelts);
end
Ints=testFaces(coeffs,T,k,formulas{3});
qxD=Ints{1};
qyD=Ints{2};
qzD=Ints{3};
ud=Ints{4};
nx=T.normals(:,[1 4 7 10]);
ny=T.normals(:,[2 5 8 11]);
nz=T.normals(:,[3 6 9 12]);
nx=nx./T.area(T.facebyele);
ny=ny./T.area(T.facebyele);
nz=nz./T.area(T.facebyele);
QxD=zeros(d2,4,Nelts);
QyD=QxD;
QzD=QxD;
UD =QxD;
Pqx=zeros(d3,Nelts);
Pqy=Pqx;
Pqz=Pqx;
Pu =Pqx;
parfor i=1:Nelts
QxD(:,:,i)=bsxfun(@times,qxD(:,T.facebyele(i,:)),nx(i,:));
QyD(:,:,i)=bsxfun(@times,qyD(:,T.facebyele(i,:)),ny(i,:));
QzD(:,:,i)=bsxfun(@times,qzD(:,T.facebyele(i,:)),nz(i,:));
UD(:,:,i)=bsxfun(@times,ud(:,T.facebyele(i,:)),tau(:,i)');
qxv=QxD(:,:,i);
qyv=QyD(:,:,i);
qzv=QzD(:,:,i);
uv=UD(:,:,i);
rhs=[qxP(:,i);
qyP(:,i);
qzP(:,i);
uP(:,i);
qxv(:)+qyv(:)+qzv(:)+uv(:)];
vect=M(:,:,i)\rhs;
Pqx(:,i)=vect(block3(1));
Pqy(:,i)=vect(block3(2));
Pqz(:,i)=vect(block3(3));
Pu(:,i)=vect(block3(4));
end
return