Skip to content
Snippets Groups Projects
Commit 863d7938 authored by graeser's avatar graeser
Browse files

Remove stabilizations

parent 90e14570
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
fd/fd_transport_stabilized.m deleted 100644 → 0
+ 0
− 70
View file @ 90e14570
function A = fd_transport_stabilized(epsilon, beta, T, n, timesteps, varargin)
tau = T/timesteps;
h = 1/(n-1);
% Peclet number
P = abs(beta*h/(2*epsilon));
if (P>1)
B = - epsilon*(1+P)*fd_laplace(n) - beta*fd_gradient(n);
else
B = - epsilon*fd_laplace(n) - beta*fd_gradient(n);
end
B = speye(n,n) + tau*B;
if length(varargin)>0
bc_type=varargin{1};
else
bc_type='dirichlet';
end
if strcmpi(bc_type,'dirichlet')
% Dirichlet boundary condition
B(1,:)=0;
B(1,1)=1;
B(end,:)=0;
B(end,end)=1;
end
if strcmpi(bc_type,'neumann')
% Neumann boundary condition
B(1,:)=0;
B(1,1)=-1/h;
B(1,2)=1/h;
B(end,:)=0;
B(end,end)=1/h;
B(end,end-1)=-1/h;
end
u0 = @(x) (x>=.4).*(x<=.5);
xx = linspace(0,1,n)';
u=xx*0;
u(:)=u0(xx);
u(1) = 0;
u(end) = 0;
plot(xx,u);
a=axis();
for k=1:timesteps
pause(.1);
rhs = u;
if strcmpi(bc_type,'neumann')
% Neumann boundary condition
rhs(1) = 0;
rhs(end) = 0;
end
u = B\rhs;
plot(xx,u);
axis(a);
end
fe/assemblers/localDiffusionConvectionAssemblerStreamlineDiffusion.m deleted 100644 → 0
+ 0
− 61
View file @ 90e14570
function A = localDiffusionConvectionAssemblerStreamlineDiffusion(grid,elementIndex,globalBasis, epsilon, beta, beta_order)
dim = size(grid.nodes,1);
order = globalBasis.order(grid, elementIndex);
n = globalBasis.localSize(grid, elementIndex);
dT = elementTransformationJacobian(grid, elementIndex);
dT_inv = inv(dT);
integrationElement = abs(det(dT));
% obtain quadrature rule of appropiate order
Q = simplexQuadratureRule(dim, max(2*order-1+beta_order, 2*(order-1)));
% evaluate basis functions at quadrature points
values = zeros(n, Q.size);
for i=1:n
values(i,:) = globalBasis.evaluate(grid, elementIndex, i, Q.points);
end
% evaluate gradients of basis functions at quadrature points
gradients = zeros(n, dim, Q.size);
for i=1:n
gradients(i,:,:) = dT_inv' * globalBasis.evaluateJacobian(grid, elementIndex, i, Q.points);
end
% evaluate beta at quadrature points
if isa(beta, 'function_handle')
% beta is either a globally defined function
values_beta = beta(elementTransformation(grid, elementIndex, Q.points, dT));
elseif (size(beta,2) == size(grid.elements))
% ... or a set a vectors given per element, i.e., a piecewise constant vector field
values_beta = repmat(beta(:,elementIndex), 1, Q.size);
else
% ... or a single vector, i.e., a constant vector field
values_beta = repmat(beta, 1, Q.size);
end
h = elementSize(grid, elementIndex);
normBeta = sum(values_beta.^2,1).^.5;
P = max(normBeta)*h/(2*epsilon);
A = zeros(n, n);
for k=1:Q.size
weight = Q.weights(k)*integrationElement;
for i=1:n
for j=1:n
% diffusion term
A(i,j) = A(i,j) + epsilon*dot(gradients(j,:,k),gradients(i,:,k)) * weight;
% convection term
A(i,j) = A(i,j) - dot(values_beta(:,k),gradients(i,:,k)) * values(j,k) * weight;
% stabilization term
normedBetaGradLambda_i = dot(values_beta(:,k),gradients(i,:,k))/normBeta(k);
normedBetaGradLambda_j = dot(values_beta(:,k),gradients(j,:,k))/normBeta(k);
A(i,j) = A(i,j) + normedBetaGradLambda_i*normedBetaGradLambda_j*epsilon*P*weight;
end
end
end
fe/fe_transport_stabilized.m deleted 100644 → 0
+ 0
− 84
View file @ 90e14570
function [u, grid] = fe_transport_stabilized(epsilon, beta, T, N, timesteps, varargin)
disp('create grid');
grid = uniformGrid(N(1),N(end));
%grid = load('circle');
refinements = 0;
for i=1:refinements
grid = refineGrid(grid);
end
disp('creating basis')
basis = lagrangeBasis(2,1);
n = basis.size(grid);
% initial value
u0 = @(x) ((x(1,:)-0.5).^2+(x(2,:)-.5).^2).^.5<.1;
tau = T/timesteps;
disp('assemble mass matrix');
MassMatrix = globalOperatorAssembler(grid, basis, @localMassAssembler, 8);
h = 1/(N-1);
P = norm(beta)*h/(2*epsilon)
disp('assemble diffusion convection matrix');
if (P>1)
DiffusionConvectionMatrix = globalOperatorAssembler(grid, basis, @(gg,ee,bb)localDiffusionConvectionAssembler(gg,ee,bb,epsilon*(1+P),beta, 0), 8);
else
DiffusionConvectionMatrix = globalOperatorAssembler(grid, basis, @(gg,ee,bb)localDiffusionConvectionAssembler(gg,ee,bb,epsilon,beta, 0), 8);
end
B = MassMatrix + tau*DiffusionConvectionMatrix;
for k=1:timesteps
pause(.1);
if (k==1)
disp('assemble right hand side for first time step');
rhs = globalFunctionalAssembler(grid, basis, @(gg,ee,bb)localL2FunctionalAssembler(gg,ee,bb,u0,1));
else
rhs = MassMatrix*u;
end
u = B\rhs;
disp(['plot solution in timepstep k=' num2str(k)]);
plotGridFunction(grid, basis, u);
if (k==1)
a=axis();
end
% axis(a);
end
%disp('incorporate boundary values');
%AA = A;
%for i=1:size(grid.nodes,2)
% if (grid.boundary(i))
% AA(i,:) = 0;
% AA(i,i) = 1;
% b(i) = 0;
% end
%end
%vtk_trimesh(grid.elements, grid.nodes(1,:)', grid.nodes(2,:)', u, u);
%vtk_trisurf(grid.elements, grid.nodes(1,:)', grid.nodes(2,:)', u, u);
fe/fe_transport_stabilized_streamline.m deleted 100644 → 0
+ 0
− 77
View file @ 90e14570
function [u, grid] = fe_transport_stabilized_streamline(epsilon, beta, T, N, timesteps, varargin)
disp('create grid');
grid = uniformGrid(N(1),N(end));
%grid = load('circle');
refinements = 0;
for i=1:refinements
grid = refineGrid(grid);
end
disp('creating basis')
basis = lagrangeBasis(2,1);
n = basis.size(grid);
% initial value
u0 = @(x) ((x(1,:)-0.5).^2+(x(2,:)-.5).^2).^.5<.1;
tau = T/timesteps;
disp('assemble mass matrix');
MassMatrix = globalOperatorAssembler(grid, basis, @localMassAssembler, 8);
disp('assemble diffusion convection matrix');
DiffusionConvectionMatrix = globalOperatorAssembler(grid, basis, @(gg,ee,bb)localDiffusionConvectionAssemblerStreamlineDiffusion(gg,ee,bb,epsilon,beta, 0), 8);
B = MassMatrix + tau*DiffusionConvectionMatrix;
for k=1:timesteps
pause(.1);
if (k==1)
disp('assemble right hand side for first time step');
rhs = globalFunctionalAssembler(grid, basis, @(gg,ee,bb)localL2FunctionalAssembler(gg,ee,bb,u0,1));
else
rhs = MassMatrix*u;
end
u = B\rhs;
disp(['plot solution in timepstep k=' num2str(k)]);
plotGridFunction(grid, basis, u);
if (k==1)
a=axis();
end
% axis(a);
end
%disp('incorporate boundary values');
%AA = A;
%for i=1:size(grid.nodes,2)
% if (grid.boundary(i))
% AA(i,:) = 0;
% AA(i,i) = 1;
% b(i) = 0;
% end
%end
%vtk_trimesh(grid.elements, grid.nodes(1,:)', grid.nodes(2,:)', u, u);
%vtk_trisurf(grid.elements, grid.nodes(1,:)', grid.nodes(2,:)', u, u);
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment