Commit 75ffd182 authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Added FEM spatial derivative matrix

parent 3c2f3852
function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
function [M, K, L, P, PHI, dphi, Dz] = femMatrices(spatialDomain, optArgs)
%FEMMATRICES computes the discretization matrices obtained by FE-Method
% [M, K, L, P, PHI, DPHI] = femMatrices(SPATIALDOMAIN, OPTARGS) computes the matrices required for
% the approximation of a system using FEM. For instance, the rhs of a pde
......@@ -8,12 +8,16 @@ function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
% The mass-matrix M is computed using
% M = int phi(z) phi(z)
% The stiffness matrix K is computed with
% K0 = int phi(z) gamma(z) phi(z) K1 = int phi(z) beta(z) (dz phi(z)) K2 = int - (dz phi(z))
% alpha(z) (dz phi(z)) K = K0 + K1 - K2
% K0 = int phi(z) gamma(z) phi(z)
% K1 = int phi(z) beta(z) (dz phi(z))
% K2 = int - (dz phi(z)) alpha(z) (dz phi(z))
% K = K0 + K1 - K2
% The input matrix L:
% L = int( phi(z) * b(z) )
% The output matrix P:
% P = int( c(z) * phi'(z) )
% Derivative Matrix
% Dz = int( phi(z) * dz phi'(z) dz
%
% PHI is the trial function used for the discretization and DPHI is the derivative of the trial
% function.
......@@ -102,12 +106,14 @@ dphidphi = reshape(dphi * dphi', 1, 4);
K0 = zeros(2,2);
K1 = zeros(2,2);
K2 = zeros(2,2);
Dz_k = zeros(2,2);
L0 = zeros(2,0);
P0 = zeros(0,2);
% initialize the assembly matrices:
M = zeros(N, N);
K = zeros(N, N);
Dz = zeros(N, N);
L = zeros(N, p);
P = zeros(m, N);
......@@ -141,7 +147,7 @@ for k = 1:N-1
if flags.alpha
% compute the elements for K2 = int dz phi(z) * alpha(z) * dz phi'(z)
% as the matrix dzphi*dzphi' is constnat the multiplication is simplified by vectorization
% as the matrix dzphi*dzphi' is constant the multiplication is simplified by vectorization
% of dzphi*dzphi' and a reshape to obtain the correct result.
k2 = alf.on(zk) * dphidphi;
K2 = reshape( numeric.trapz_fast(zk, k2'), 2, 2);
......@@ -159,12 +165,19 @@ for k = 1:N-1
P0 = reshape( numeric.trapz_fast(zk, c0'), m, 2);
end
if nargout >= 7
Dz_ = misc.multArray(phi, repmat( dphi', [optArgs.Ne,1]), 3, 3, 1);
Dz_k = reshape( numeric.trapz_fast(zk, reshape(Dz_, element.n, 4, 1)' ), 2, 2);
end
% assemble the matrices with the elements:
M(idxk, idxk) = M(idxk, idxk) + M0 * zDiff(k);
K(idxk, idxk) = K(idxk, idxk) + K0 + K1 - K2;
L(idxk, :) = L(idxk, :) + L0;
P(:, idxk) = P(:, idxk) + P0;
Dz(idxk, idxk) = Dz(idxk, idxk) + Dz_k;
end
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment