Commit fbce5bec authored by Ferdinand Fischer's avatar Ferdinand Fischer
Browse files

Correcting the documentation regarding the new input parameters.

parent 0e13619e
function [M, K, L, P, PHI, dphi] = femMatrices(spatialDomain, optArgs)
%FEMMATRICES computes the discretization matrices obtained by FE-Method
% [M, K, L, P, PHI] = femMatrices(varargin) computes the matrices required
% for the approximation of a system using FEM. For instance, the rhs of a pde
% dz (alpha(z) dz x(z)) + beta(z) dz x(z) + gamma(z) x(z) + b(z) u
% y = c(z) x(z)
% is approximated as K x(z)
% The matrices that will be returned as described in the following:
% [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
% PDE[ x(z) ] = dz (alpha(z) dz x(z)) + beta(z) dz x(z) + gamma(z) x(z) + b(z) u
% y = c(z) x(z)
% is approximated as K x(z) The matrices that will be returned are:
% 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) )
%
% PHI is the trial function used for the discretization.
% PHI is the trial function used for the discretization and DPHI is the derivative of the trial
% function.
%
% The argument SPATIALDOMAIN is mandatory and must be a quantity.Domain object. The further system
% parameters as
% "alpha", "beta", "gamma", "b", "c"
% must be quantity.Discrete objects. With the parameter
% "Ne"
% as positive integer number one can specify the number of discretization points for one element.
% The parameter
% "silent"
% is a logical flag to supress the progressbar.
%
% The input arguments can be given as name-value-pairs. The parameters
% are described in the following, the default values are written in (.).
% 'grid' (linspace(0, 1, 11)): grid of resulting finite elements
% 'Ne'(111): the implementation uses numerical integration for each element.
% The parameter 'Ne' defines the grid for this integration.
% 'alpha' : the system parameter alpha, should be a quantity.Discrete.
% 'beta' : the system parameter beta, should be a quantity.Discrete
% 'gamma' : the system parameter gamma, should be a quantity.Discrete
% 'b' : the input vector, should be a quantity.Discrete
% 'c' : the output vector, should be a quantity.Discrete
arguments
spatialDomain (1,1) quantity.Domain;
......@@ -42,7 +40,7 @@ arguments
optArgs.silent logical = false;
end
% set the system properties
% set the system properties to local variables, for better readability and (maybe) faster access
z = spatialDomain.grid;
N = spatialDomain.n;
z0 = spatialDomain.lower;
......@@ -65,9 +63,9 @@ element = quantity.EquidistantDomain("z", 0, (z1 - z0) / (N - 1), ...
'stepNumber', optArgs.Ne);
%% verify that grid is homogenious
% Yet, this function only supports homogenious grids for most output parameters.
% Only, the mass matrix M is calculated correct for nonhomogenious grids. The
% following line ensures, that if the grid is not homogenious, only M is returned.
% Yet, this function only supports homogenious grids for most output parameters. Only, the mass
% matrix M is calculated correct for nonhomogenious grids. The following line ensures, that if the
% grid is not homogenious, only M is returned.
isGridHomogenious = all(diff(z) - element.upper < 10*eps);
if ~isGridHomogenious
assert(nargout <= 1, 'non-homogenious grids only supported for mass-matrix');
......@@ -92,8 +90,8 @@ zDiff = diff(z);
%% compute the stiffness matrix
% initialize the trial functions for the computation for element matrices.
% for further comments, see the description in the loop
% initialize the trial functions for the computation for element matrices. for further comments, see
% the description in the loop
pp = reshape(double(PHI * PHI'), element.n, 4);
phi = double(PHI);
dphi = [ -1 / element.upper; ...
......@@ -127,27 +125,24 @@ for k = 1:N-1
if flags.gamma
% compute the elements for K0 = int phi(z) * gamma(z) * phi'(z)
% the matrix phi(z) * phi'(z) is the same for each element. Thus,
% this matrix is computed a priori. To simplify the multiplication
% with the scalar value gamma(z), the matrix phi*phi' is vectorized
% and the result of phi(z)phi'(z) alpha(z) is reshaped as a matrix.
% the matrix phi(z) * phi'(z) is the same for each element. Thus, this matrix is computed
% a priori. To simplify the multiplication with the scalar value gamma(z), the matrix
% phi*phi' is vectorized and the result of phi(z)phi'(z) alpha(z) is reshaped as a matrix.
k0 = misc.multArray(pp, gma.on(zk), 3, 2, 1);
K0 = reshape( numeric.trapz_fast(zk, k0'), 2, 2);
end
if flags.beta
% compute the elements for K1 = int phi(z) * beta(z) * dz phi'(z))
% as the values dz phi(z) are constants, the multiplication can be.
% simplified.
% as the values dz phi(z) are constants, the multiplication can be. simplified.
k1 = misc.multArray(phi, bta.on(zk) * dphi', 3, 3, 1);
K1 = reshape( numeric.trapz_fast(zk, reshape(k1, element.n, 4, 1)' ), 2, 2);
end
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 of dzphi*dzphi' and a reshape to
% obtain the correct result.
% as the matrix dzphi*dzphi' is constnat 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);
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