Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (master) Browse

fennecs/matlab/splinepressure.m
dd5453d97d74master

Last Change
Raw File

splinepressure.m
View Options

classdef splinepressure < h5quantity
properties
invVolume
knotsr
knotsz
femorder
postscale
end
methods (Access=public)
function obj=splinepressure(filename, dataset, knotsr, knotsz, femorder, Volume, vnorm, postscale)
if nargin<7
vnorm=1;
end
obj=obj@h5quantity(filename, dataset, length(knotsr)-2*femorder(2), length(knotsz)-2*femorder(1), vnorm);
obj.knotsr=knotsr;
obj.knotsz=knotsz;
obj.femorder=double(femorder);
obj.invVolume=1./Volume;
obj.invVolume(isinf(obj.invVolume))=0;
if nargin <8
obj.postscale=1;
else
obj.postscale=postscale;
end
end
function quantity=coeffs(obj,indices)
if strcmp(indices{1},':')
ij=1:6;
else
ij=indices{1};
end
if strcmp(indices{2},':')
r=1:obj.nr+obj.femorder(2)-1;
else
r=indices{2};
end
if strcmp(indices{3},':')
z=1:obj.nz+obj.femorder(1)-1;
else
z=indices{3};
end
if strcmp(indices{4},':')
t=1:obj.nt;
else
t=indices{4};
end
quantity=zeros(length(ij),(obj.nr+obj.femorder(2)-1),(obj.nz+obj.femorder(1)-1),length(t));
%temp=zeros((obj.nz+obj.femorder(1)-1)*(obj.nr+obj.femorder(2)-1),length(t));
for tempi=1:length(ij)
temp=obj.readtensorindex(ij(tempi),t);
for timei=1:size(temp,2)
quantity(tempi,:,:,timei)=permute(reshape(temp(:,timei),obj.nz+obj.femorder(1)-1,obj.nr+obj.femorder(2)-1),[2,1,3]).*obj.invVolume;
end
end
quantity=quantity*obj.scale;
end
function quantity=val(obj,indices)
if strcmp(indices{1},':')
ij=1:6;
else
ij=indices{1};
end
if strcmp(indices{2},':')
r=1:obj.nr;
else
r=indices{2};
end
if strcmp(indices{3},':')
z=1:obj.nz;
else
z=indices{3};
end
if strcmp(indices{4},':')
t=1:obj.nt;
else
t=indices{4};
end
count=length(t);
temp=obj.coeffs({ij,':',':',t});
quantity=zeros(length(ij),max(r)-min(r)+1,max(z)-min(z)+1,count);
for j=1:size(temp,1)
for i=1:size(temp,4)
quantity(j,:,:,i)=obj.postscale.*fnval(spmak({obj.knotsr,obj.knotsz},squeeze(temp(j,:,:,i))),{obj.knotsr(r+obj.femorder(2)),obj.knotsz(z+obj.femorder(1))});
end
end
end
function quantity=der(obj,indices)
if strcmp(indices{1},':')
ij=1:6;
else
ij=indices{1};
end
if strcmp(indices{2},':')
r=1:obj.nr;
else
r=indices{2};
end
if strcmp(indices{3},':')
z=1:obj.nz;
else
z=indices{3};
end
if strcmp(indices{4},':')
t=1:obj.nt;
else
t=indices{4};
end
if length(indices)<5
order=[1,1];
else
order=indices{5};
end
count=length(t);
temp=obj.coeffs({ij,':',':',t});
quantity=zeros(length(ij),max(r)-min(r)+1,max(z)-min(z)+1,count);
for j=1:size(temp,1)
for i=1:size(temp,4)
quantity(j,:,:,i)=fnval(fnder(...
spmak({obj.knotsr,obj.knotsz},squeeze(temp(j,:,:,i)))...
,order),{obj.knotsr(r+obj.femorder(2)),obj.knotsz(z+obj.femorder(1))});
end
end
end
function ind=end(obj,k,n)
switch k
case 1
ind=6;
case 2
ind=obj.nr;
case 3
ind=obj.nz;
case 4
ind=obj.nt;
case default
error('Invalid number of dimensions');
end
end
end
methods (Access=private)
function quantity = readtensorindex(obj, ij, t)
if strcmp(ij,':') || length(ij)>1
error('Unable to read several indices at the same time');
end
if strcmp(t,':')
t=1:obj.nt;
end
ui=zeros((obj.nz+obj.femorder(1)-1)*(obj.nr+obj.femorder(2)-1),length(t));
uj=ui;
n=ui;
vij=ui;
switch ij
case 1
i=1;
j=1;
case 2
i=1;
j=2;
case 3
i=1;
j=3;
case 4
i=2;
j=2;
case 5
i=2;
j=3;
case 6
i=3;
j=3;
case default
error('Invalid Pressure index');
end
for k=1:length(t)
vij(:,k) = h5read(obj.filename, obj.dataset,[ij+4 1 t(k)],[1 Inf 1]);
ui(:,k) = h5read(obj.filename, obj.dataset,[i+1 1 t(k)],[1 Inf 1]);
uj(:,k) = h5read(obj.filename, obj.dataset,[j+1 1 t(k)],[1 Inf 1]);
n(:,k) = h5read(obj.filename, obj.dataset,[1 1 t(k)],[1 Inf 1]);
end
n=1./n;
n(isinf(n))=0;
quantity=vij-ui.*uj.*n;
end
end
end
c4science · Help