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

Add streamline stabilization

parent 84a99a28
Branches
No related tags found
No related merge requests found
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
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% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment