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

Initial commit

parents
Branches
No related tags found
No related merge requests found
Showing
with 586 additions and 0 deletions
function y = globalFunctionalAssembler(grid, basis, localFunctionalAssembler)
% dimension of ansatz space
n = basis.size(grid);
% allocate vector
y = zeros(n,1);
% loop over all elements
for e=1:size(grid.elements,2)
% compute local stiffness matrix
localB = localFunctionalAssembler(grid,e,basis);
% compute global indices
indices = basis.index(grid, e);
% add local vector to corresponding global vector entries
y(indices) = y(indices) + localB;
end
function matrix = globalOperatorAssembler(grid, basis, localAssembler, average_nnz)
% dimension of ansatz space
n = basis.size(grid);
% allocate matrix
matrix= spalloc(n,n,average_nnz*n);
% loop over all elements
for e=1:size(grid.elements,2)
% compute local stiffness matrix
localMatrix = localAssembler(grid, e, basis);
% compute global indices
indices = basis.index(grid, e);
% add local matrix to corresponding global matrix entries
matrix(indices,indices) = matrix(indices,indices) + localMatrix;
end
function A = localDiffusionConvectionAssembler(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
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;
end
end
end
function A = localGradientAssembler(grid,elementIndex,globalBasis, 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, 2*order-1+beta_order);
% 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
A = zeros(n, n);
for k=1:Q.size
weight = Q.weights(k)*integrationElement;
for i=1:n
z = dot(values_beta(:,k),gradients(i,:,k))*weight;
for j=1:n
A(i,j) = A(i,j) - values(j,k) * z;
end
end
end
function y = localL2FunctionalAssembler(grid, elementIndex, globalBasis, f, quadOrder)
dim = size(grid.nodes,1);
order = globalBasis.order(grid, elementIndex);
n = globalBasis.localSize(grid, elementIndex);
dT = elementTransformationJacobian(grid, elementIndex);
integrationElement = abs(det(dT));
% obtain quadrature rule of appropiate order
Q = simplexQuadratureRule(dim, order+quadOrder);
% 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 f at quadrature points
values_f = f(elementTransformation(grid, elementIndex, Q.points, dT));
y = zeros(n,1);
for k =1:Q.size
z = values_f(k)*Q.weights(k)*integrationElement;
for i=1:n
y(i) = y(i) + values(i,k)*z;
end
end
end
function A = localLaplaceAssembler(grid,elementIndex,globalBasis)
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, (order-1)*2);
% 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
A = zeros(n, n);
for k=1:Q.size
weight = Q.weights(k)*integrationElement;
% only compute lower triangle and diagonal entries
for i=1:n
z = gradients(i,:,k)*weight;
for j=1:i
A(i,j) = A(i,j) + dot(gradients(j,:,k),z);
end
end
end
% copy upper triangle (exploit symmetrie)
for i=1:n
for j=(i+1):n
A(i,j) = A(j,i);
end
end
function A = localMassAssembler(grid,elementIndex,globalBasis)
dim = size(grid.nodes,1);
order = globalBasis.order(grid, elementIndex);
n = globalBasis.localSize(grid, elementIndex);
dT = elementTransformationJacobian(grid, elementIndex);
integrationElement = abs(det(dT));
% obtain quadrature rule of appropiate order
Q = simplexQuadratureRule(dim, order*2);
% 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
A = zeros(n, n);
for k=1:Q.size
weight = Q.weights(k)*integrationElement;
% only compute lower triangle and diagonal entries
for i=1:n
z = values(i,k)*weight;
for j=1:i
A(i,j) = A(i,j) + values(j,k)*z;
end
end
end
% copy upper triangle (exploit symmetrie)
for i=1:n
for j=(i+1):n
A(i,j) = A(j,i);
end
end
function B = lagrangeBasis(dim, order)
if ((order==1) & (dim==1))
B = struct( ...
'dim', 1, ...
'evaluate', @D1P1Evaluate, ...
'evaluateJacobian', @D1P1EvaluateJacobian, ...
'localSize', @(grid, elementIndex)(2), ...
'localInterpolation', @P1LocalInterpolation, ...
'order', @(grid, elementIndex)(1),...
'index', @P1Index,...
'size', @(grid)(size(grid.nodes,2)));
elseif ((order==1) & (dim==2))
B = struct( ...
'dim', 2, ...
'evaluate', @D2P1Evaluate, ...
'evaluateJacobian', @D2P1EvaluateJacobian, ...
'localSize', @(grid, elementIndex)(3), ...
'localInterpolation', @P1LocalInterpolation, ...
'order', @(grid, elementIndex)(1),...
'index', @P1Index,...
'size', @(grid)(size(grid.nodes,2)));
end
end
function u = P1LocalInterpolation(grid, elementIndex, f, df, ddf, dddf)
element = grid.nodes(grid.elements(:,elementIndex));
n = size(element,2);
u = zeros(1,n);
for i=1:n
u(i) = f(element(:,i));
end
end
function globalIndices = P1Index(grid, elementIndex, localIndices)
% Compute global indices of basis functions with given localIndices in element.
dim = size(grid.nodes,1);
globalIndices = grid.elements(1:(dim+1), elementIndex);
if (nargin==3)
globalIndices = globalIndices(localIndices);
end
end
function y = D1P1Evaluate(grid, elementIndex, localIndex, x)
switch localIndex
case 1,
y = x;
case 2,
y = 1-x;
end
end
function Dx = D1P1EvaluateJacobian(grid, elementIndex, localIndex, x)
Dx = x*0;
switch localIndex
case 1,
Dx(1,:) = 1;
case 2 ,
Dx(1,:) = -1;
end
end
function y = D2P1Evaluate(grid, elementIndex, localIndex, x)
switch localIndex
case 1,
y = 1 - x(1,:) - x(2,:);
case 2 ,
y = x(1,:);
case 3,
y = x(2,:);
end
end
function Dx = D2P1EvaluateJacobian(grid, elementIndex, localIndex, x)
Dx = x*0;
switch localIndex
case 1,
Dx(1,:) = -1;
Dx(2,:) = -1;
case 2 ,
Dx(1,:) = 1;
Dx(2,:) = 0;
case 3,
Dx(1,:) = 0;
Dx(2,:) = 1;
end
end
function du = gradientValues(grid, basis, u)
dim = size(grid.nodes,1);
% center of reference simplex
center = repmat(1/(dim+1),dim,1);
% allocate matrix
du= zeros(dim, size(grid.elements,2));
% loop over all elements
for e=1:size(grid.elements,2)
% compute geometry information
simplex = grid.nodes(:,grid.elements(1:(dim+1),e));
dT = zeros(dim);
for i=1:dim
dT(:,i) = simplex(:,i+1) - simplex(:,1);
end
integrationElement = abs(det(dT));
dTInvTrans = inv(dT');
% compute global indices
indices = basis.index(grid, e);
du(:,e) = 0;
for i=1:globalBasis.localSize(grid, e)
gradient_i = dTInvTrans * globalBasis.evaluateJacobian(grid, e, i, center);
du(:,e) = du(:,e) + u(indices(i))* gradient_i;
end
end
function u = interpolateFunction(grid, basis, varargin)
% dimension of ansatz space
n = basis.size(grid);
% initialize vector
u = zeros(n, 1);
% loop over all elements
for e=1:size(grid.elements,2)
E = grid.nodes(1, grid.elements(:,e));
% transform function and derivatives to element to
fLocal = {};
for i=1:length(varargin)
fLocal{i} = @(x)varargin{i}(E(1) + (E(2)-E(1))*x);
end
% compute local interpolation
uLocal = basis.localInterpolation(grid, e, fLocal{:});
% local size
nLocal = basis.localSize(grid, e);
% compute global indices
for j=1:nLocal
u(basis.index(grid, e, j)) = uLocal(j);
end
end
function plotGridFunction(grid, basis, u)
dim = size(grid.nodes,1);
if (dim==2)
% referenceElement = generateGrid([0 0; 1 0; 0 1]', [1 2 3]');
n_nodes = size(grid.nodes, 2);
trisurf(grid.elements(1:3,:)', grid.nodes(1,:)', grid.nodes(2,:)', u(1:n_nodes), u(1:n_nodes));
end
function [xx,yy] = sampleGridFunction(grid, basis, u, localPoints)
dim = size(grid.nodes,1);
% store all sampled values and produce one plot only
xx = zeros(length(localPoints), size(grid.elements,2));
yy = zeros(length(localPoints), size(grid.elements,2));
for e=1:size(grid.elements,2)
vertices = grid.elements(1:(dim+1),e);
element = grid.nodes(vertices);
xx(:,e) = linspace(element(1), element(2), length(localPoints));
for i=1:basis.localSize(grid, e)
globalIndex = basis.index(grid, e, i);
values = basis.evaluate(grid, e, i, localPoints);
yy(:,e) = yy(:,e) + (u(globalIndex) * values(:));
end
end
function T = assemble_transfer_operator(coarse_grid, fine_grid)
n_coarse_nodes = size(coarse_grid.nodes, 2);
n_fine_nodes = size(fine_grid.nodes, 2);
n_fine_elements = size(fine_grid.elements, 2);
T = spalloc(n_fine_nodes, n_coarse_nodes, 3*n_coarse_nodes);
for k = 1:n_fine_elements
end
File added
function geo = elementGeometry(grid, elementIndex)
dim = size(grid.nodes,1);
vertices = grid.nodes(:,grid.elements(1:(dim+1),elementIndex));
dT = zeros(dim);
for i=1:dim
dT(:,i) = vertices(:,i+1) - vertices(:,1);
end
geo = struct( ...
'dim', dim, ...
'vertices', vertices, ...
'jacobian', dT,
'inverseJacobian', inv(dT),
'global', @(x)D1P1EvaluateJacobian, ...
'localSize', @(grid, elementIndex)(2), ...
'localInterpolation', @P1LocalInterpolation, ...
integrationElement = abs(det(dT));
dTInvTrans = inv(dT');
function h = elementSize(grid, elementIndex)
dim = size(grid.nodes,1);
simplex = grid.nodes(:,grid.elements(1:(dim+1),elementIndex));
h = 0;
for i=1:(dim+1)
for j=(i+1):(dim+1)
h = max(h, norm(simplex(:,i)-simplex(:,j)));
end
end
function xGlobal = elementTransformation(grid, elementIndex, xLocal, dT)
firstVertex = grid.nodes(:,grid.elements(1,elementIndex));
xGlobal = xLocal*0;
for i=1:size(xLocal,2)
xGlobal(:,i) = firstVertex + dT*xLocal(:,i);
end
function xLocal = elementTransformationInverse(grid, elementIndex, xGlobal, dT_inv)
firstVertex = grid.nodes(:,grid.elements(1,elementIndex));
xLocal = xGlobal*0;
for i=1:size(xLocal,2)
xLocal(:,i) = dT_inv*(xGlobal(:,i) - firstVertex);
end
function dT = elementTransformationJacobian(grid, elementIndex)
dim = size(grid.nodes,1);
simplex = grid.nodes(:,grid.elements(1:(dim+1),elementIndex));
dT = zeros(dim,dim);
for i=1:dim
dT(:,i) = simplex(:,i+1) - simplex(:,1);
end
function grid = generateGrid(raw_nodes, raw_elements)
% Generate a grid containing edge and boundary information from raw node and element lists
% Indices are ordered such that edge i is opposite to node i in each element
%
% Usage:
% grid = generate_grid(nodes, raw_elements)
%
% Input:
% raw_nodes - a 2xn list of nodes
% raw_elements - a 3xm list of elements
%
% Output:
% grid - a structure containing:
% grid.nodes list of node given by coordinates (1:2,:) and is boundary information (3,:)
% grid.boundary is boundary information per node (1,:)
% grid.edges list of edges given by node numbers (1:2,:)
% grid.elements list of triangles given by node numbers (1:3,:) and edge numbers (4:6,:)
n_nodes = size(raw_nodes, 2);
n_elements = size(raw_elements, 2);
edge_mat = spalloc(n_nodes, n_nodes, 10*n_nodes);
if (size(raw_elements,1)==4)
elements = [raw_elements(1:3,:) ; zeros(3, n_elements); raw_elements(4,:)];
else
elements = [raw_elements(1:3,:) ; zeros(4, n_elements)];
end
%elements = [raw_elements(1:3,:) ; zeros(3, n_elements)];
n_edges = 0;
for k = 1:n_elements
t = sort(elements(1:3, k));
ij = 1;
for i=1:3
for j=(i+1):3
if (edge_mat(t(i),t(j))==0)
n_edges = n_edges + 1;
edge_mat(t(i),t(j)) = n_edges;
end
elements(3+ij, k) = edge_mat(t(i), t(j));
ij = ij+1;
end
end
end
[i,j,k] = find(edge_mat);
edges = zeros(3, n_edges);
edges(1:2, k) = [i' ; j'];
edges(3,:) = 2;
for i=1:n_elements
indices = elements(4:6,i);
edges(3, indices) = edges(3,indices)-1;
elements(1:6,i) = sort_indices(elements(1:6,i), edges);
end
boundary = zeros(1,size(raw_nodes,2));
boundary_edges = find(edges(3,:)==1);
boundary(edges(1,boundary_edges)) = 1;
boundary(edges(2,boundary_edges)) = 1;
%raw_nodes(3,:) = boundary;
grid = struct('nodes', raw_nodes, 'elements', elements, 'boundary', boundary, 'edges', edges);
end
function t = sort_indices(t,edges)
t = [sort(t(1:3)); t(4:6)];
for i=1:2
for j=(i+1):3
if ((edges(1,t(3+i))==t(i)) || (edges(2,t(3+i))==t(i)))
dummy = t(3+i);
t(3+i) = t(3+j);
t(3+j) = dummy;
end
end
end
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment