% date: 16.12.04
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simple FEM program
%
% schreiben Sie in die Funktionen
%
% fems, femm, femc, feme und femf
%
% die fehlenden Anweisungen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = fem
%----------------
% compute FEM solution U of the problem
%--------------------------------------
global TEST;
%show results (1=show, 0=don't)
global ourTEST
ourTEST = [0 0 0 0 0]; %femf, feme, femc, femm, fems
%TEST = 1; % test display
TEST = 0; % no test display
% define/ init problem data
%bnr = 4;
%ng = 1;
%coeffall(bnr,ng); % init example
femi; % initialize fem computation
femc; % compile SV and SM
femm; % modify SV
fems; % compute solution U
femp; % plot solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% functions for FEM computation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [] = fems
%-----------------
% solve system equation SM * U = SV
%----------------------------------------------------
global L N
global Z
global SV SM
global U % calculated
global ourTEST
global TEST
global ERROR
global maxERROR
%solve equation
U = SM\SV;
U = reshape(U,[L N]);
ERROR = [NaN]; maxERROR = NaN; exactResult = NaN;
%calc Error if L=1
if L==1
% calc exact result
for i = 1:N
V(i) = uexact([Z(i,1) Z(i,2)]);
end
exactResult = reshape(V,[1 N]);
% calc ERROR
for i = 1:N
ERROR( i ) = abs( V(i) - U(i) );
end
ERROR = reshape( ERROR,[1 N] );
maxERROR = max( ERROR );
end
%Show Results
if ourTEST(5)
disp('fems ready, U and ERROR calculated:');
disp('U = ');disp(U);
disp('exact result ');disp(exactResult);
disp('ERROR = '); disp(ERROR);
disp('maxERROR = '); disp(maxERROR);
disp('SYSTEM EQUATION SM * U = SV is solved for U');
end
function [] = femm
%-----------------
% MODIFY SV and SM for DIRICHLET boundary condition
%
% use variable ID = 1 x N - matrix of kinnd of knotes
%----------------------------------------------------
global SV % modified
global SM % modified
global ID % Nodetable
global Z % x,y values
global L
global ourTEST
global TEST
ind = find(ID == 1); % all dirichlet nodes k
for k = ind %dirichlet node k
uk = u0(Z(k,:));
indk = ((L*(k-1)) + 1) : (L*(k-1)+L); % Columnnumber
% overwrite SV with new and old known values
SV = SV - SM(:,indk) * uk;
SV(indk) = uk;
% modify SM
SM(indk,:) = 0;
SM(:,indk) = 0;
SM(indk,indk) = eye(L);
end
figure; spy(SM); figure
%Debug
if ourTEST(4)
disp('femm ready, SV modified: ');
disp(SV);
disp('--------------------- press any key --------------------');
pause
end
function [] = femc
%-----------------
% compile global system variables
%
% SM = C = L x L cell = stiffness matrices for each solution component
% SV = b = L x 1 cell = load vectors
%
% some other used variables:
%
% L = # of unknown functions = components of u
% M = # of domain triangles
% TOP = all domain knotes
%
% SME = system matrix on domain ELEMENT # r
% SVE = system vector on domain ELEMENT # r
%
% call function feme(r) = compute elemnt variables SME and SVE
%--------------------------------------------------------------
global TOP
global L M N
global SVE SME
global SV SM % calculated
global ourTEST
SV = zeros(N*L,1); %b:=0
SM = zeros(N*L,N*L); %C:=0
for r =1:M
feme(r); %calc SVE and SME
for alpha = 1:3
i = TOP(r,alpha);
indi = L*(i-1)+1 : L*(i-1)+L;
indalpha = L*(alpha-1)+1 : L*(alpha-1)+L;
SV(indi) = SV(indi) + SVE(indalpha);
for beta = 1:3
k = TOP(r,beta);
indk = L*(k-1)+1 : L*(k-1)+L;
indbeta = L*(beta-1)+1: L*(beta-1)+L;
SM(indi,indk) = SM(indi,indk) + SME(indalpha,indbeta);
end
end
end
figure; spy(SM);
if ourTEST(3)
disp('femc ready, SV and SM calculated:');
disp('SV:');disp(reshape(SV, [L N]));
disp('press any key');pause;
disp('SM:');disp(SM);
disp('--------------------- press any key --------------------');
pause;
end
function [] = feme(r)
%--------------------
% compute on DOMAIN triangle/element # r the system variables
%
% SVE = b = L x 1 cell = lokal ELEMENT load vector
% SME = C = L X L cell = lokal ELEMENT stiffness matrix
%
% where
%
% r = domain # of actual triangle/element
%
% use FFIA/FFIB = cells of FF integrals for element r
%
% FFIA = {FIA1 FIA2 FIAx FIAy FIAxx FIAxy FIAyx FIAyy} = area integrals
% FFIB = {FIB1 FIB2 FIB3} = boundary integrals
%
% call function femf(r) = compute FF integrals FFIA and FFIB on element r
%----------------------------------------------------------------------------
global ourTEST
global TEST
global INC
global L
global XE
global FFIA FFIB
global SVE SME % calculated
femf(r); % compute FF integrals FFIA and FFIB on element r
% init SVE,SME
SVE=zeros( 3*L,1 );
SME=zeros( 3*L,3*L );
% get the center of element r
center= zeros( 1,2 );
for i=1:3
center=center+XE( i,: );
end
center=center/3;
% get Coeffs for AreaInt1&2
F1=f1( center );
F2=f2( center );
SVE = kron( FFIA{1},f0( center ) );
SME = kron( FFIA{2}, a0(center) );
SME = SME + kron( FFIA{3}, F1{1} );
SME = SME + kron( FFIA{4}, F2{1} );
SME = SME + kron( FFIA{5}, F1{2} );
SME = SME + kron( FFIA{7}, F2{2} );
SME = SME + kron( FFIA{6}, F1{3} );
SME = SME + kron( FFIA{8}, F2{3} );
if TEST && ourTEST(2)
disp('r = ');disp(r);
disp('SVE:');disp(SVE);
disp('SME:');disp(SME);
end
for delta=1:3
if INC( r,delta )>0
% compute barycenter
barycenter = ( 3*center-XE( delta,: ) ) / 2;
% determine which border and the corresponding function g and a
if INC( r,delta )==1
my_g=g1( barycenter );
my_alpha=zeros( L,L );
else
my_g=g2( barycenter );
my_alpha=a2( barycenter );
end
% adjust load vector and stiffness matrix
SVE = SVE + kron( sum( FFIB{delta},2 ), my_g );
SME = SME + kron( FFIB{delta}, my_alpha );
if TEST && ourTEST(2)
disp('delta-Iteration');disp(delta);
disp('SVE:');disp(SVE);
disp('SME:');disp(SME);
end
end
end
if ourTEST(2)
disp('feme ready, SVE and SME calculated:');
disp('r =');disp(r);
disp('SVE:');disp(SVE);
disp('SME:');disp(SME);
disp('--------------------- press any key --------------------');
pause;
end
function [] = femf(r)
%--------------------
% calculate element dependent integrals
% FFIA = AREA and FFIB = BOUNDARY INTEGRALS
% of products of FF on triangle element number r
%
% FFIA = {FIA1 FIA2 FIAx FIAy FIAxx FIAxy FIAyx FIAyy} = cell of AREA integrals
% FIA1,FIA2 = integrals of 1 and 2 FF with no derivation
% FIAx,FIAy = integrals of 2 FF with one derivation
% FIAxx,FIAxy,FIAyx,FIAyy = integrals of 2 FF with two derivations
%
% FFIB = {FIB1 FIB2 FIB3} = cell of BOUNDARY integrals
%
% some used/defined variables:
%
% JD = jacobi determinante JA = adjunct of jacobi matrix
% XE = element knotes = 3 x 2 - matrix
%
% FFSIA/FFSIB = STANDARD area/boundary integrals
%------------------------------------------------------------
global ourTEST
global TEST
global Z TOP
global FFSIA FFSIB % given
global XE
global FFIA FFIB % calculated
if TEST && ourTEST(1)
disp( 'r= ' );
disp( r );
end
% element nodes
XE = zeros( 3,2 );
XE = [ Z( TOP( r,1 ),: ); Z( TOP( r,2 ),: ); Z( TOP( r,3 ),: )];
% calculate the jacobi matrix
JM = transpose( [ XE(2,:)-XE(1,:); XE(3,:)-XE(1,:)] );
%debug jacobi matrix
if TEST && ourTEST(1)
disp( 'This is JM' );
disp( JM );
end
% calculate JD
JD = det( JM );
% calculate JA
JA = [ JM( 2,2 ) -JM( 1,2 ); -JM( 2,1 ) JM( 1,1 ) ];
% element lengths
XE_diff = [ XE(3,:)-XE(2,:);
XE(1,:)-XE(3,:);
XE(2,:)-XE(1,:)];
% calculate the Ls (don't use "L" as Var because needed elsewhere -> "eL")
eL = [ norm( XE_diff( 1,: ) ) norm( XE_diff( 2,: ) ) norm( XE_diff( 3,: ) ) ];
%calcuöate FIA1 FIA2
FFIA{1} = JD * FFSIA{1};
FFIA{2} = JD * FFSIA{2};
% Area-Integrals
n=1;
for i=1:2
FFIA{i+2} = JA(1,i) * FFSIA{3} + JA(i,2) * FFSIA{4};
for k=1:2
% calculate Sum
FFIA{n+4} = 0;
m=1;
for s=1:2
for t=1:2
FFIA{n+4} = FFIA{n+4} + JA(s,i) * JA(t,k) * FFSIA{m+4};
m = m+1;
end
end
FFIA{n+4} = FFIA{n+4}/JD;
n = n+1;
end
end
% Border-Integrals
for delta=1:3
FFIB{delta} = eL(delta)*FFSIB{delta};
end
if ourTEST(1)
disp('femf ready, FFIA and FFIB calculated:');
disp('r =');disp(r);
disp('FFIA:');disp(FFIA);
disp('FFIB:');disp(FFIB);
disp('--------------------- press any key --------------------');
pause;
end
function [] = femi
%-----------------
% initialize variables needed for fem computation; compute
%
% M = # of domain triangles
% N = # of all knotes
%
% FFSIA/FFSIB = STANDARD area/boundary integrals
% FFSIA = {I1 I2 Ix Iy Ixx Ixy Iyx Iyy}
% 1 2 3 4 5 6 7 8
% FFSIB = {I1 I2 I3}
%
% ID = 1 x N - matrix of kinnd of knotes
% INC = M x 3 - matrix of kiind of triangle sides
%
% call function setbart = define ID and INC
% ploted = plot triangles
%----------------------------------------------------------
global Z TOP
global L M N
global FFSIA FFSIB % calculated
global ID
[M dummy] = size(TOP); % M = # of triangles
[N dummy] = size(Z); % N = # of global knotes
disp(['[L M N] = ' num2str([L M N])])
% plot domain tringles
%---------------------
figure(1); clf; % clear figure 1 for plot
if M < 100 % M = # of domain triangles
plotel(Z,[1 2]); % plot triangles with triangle #
else
plotel(Z,[0 1]); % plot triangles without number
end
disp('domain triangles are ploted. enter any key'); pause
% define standard area and boundary integrals
I1 = [1; 1; 1] / 6;
I2 = [2 1 1; 1 2 1; 1 1 2]/24;
Ix = [-1 -1 -1; 1 1 1; 0 0 0]/6;
Iy = [-1 -1 -1; 0 0 0; 1 1 1]/6;
Ixx = [1 -1 0; -1 1 0; 0 0 0]/2;
Ixy = [1 0 -1; -1 0 1; 0 0 0]/2; Iyx = Ixy';
Iyy = [1 0 -1; 0 0 0; -1 0 1]/2;
FFSIA = {I1 I2 Ix Iy Ixx Ixy Iyx Iyy}; % standard area integrals
%I = [0 1 1; 1 0 1; 1 1 0]/2;
I1 = [0 0 0; 0 2 1; 0 1 2]/6;
I2 = [2 0 1; 0 0 0; 1 0 2]/6;
I3 = [2 1 0; 1 2 0; 0 0 0]/6;
FFSIB = {I1 I2 I3}; % standard boubdary integrals
setbart; % calculate INC = M x 3 matrix of kind of boundary conditions of all triangle sides
% ID = 1 x N matrix of kind of all global knotes
%ID(7) = 0; % set some extra dirichlet knotes
%ind = find(ID == 1)
%disp(['all dirichlet knotes']); pause
function [] = setbart
%--------------------
% set/calculate/initialize ID and INC
% ID = 1 x N - matrix of DIRICHLET knotes
% ID(i) = 1 knot i is a DIRICHLET knot
% = -2 knot i is inner knote
% = 0 knot i is boundary knot and not dirichlet
% INC = M x 3 matrix
% = -2 for all INNER sides and
% = art of BND for all BOUNDARY sides
%
% INC(r,jj) = bart <=> side jj in element r is of type bart
%
% side jj = 1 = (2 3) , jj = 2 = (3 1) and jj = 3 = (1 2)
%
% bart > -2 <=> jj = BOUNDARY side for component l
% = -1 = homogen NC 0 = Dirichlet 1 = Neumann 2 = Cauchy
%
% used global variable:
%
% BART = given matrix of boundary types (defied in example files)
%-----------------------------------------------------------------------
global TOP M N
global BART % given boundary types
global INC ID % are calculated !!!!!!!!
global TEST
INC = (-2)*ones(M,3); % initialize INC = -2 for all triangle sides
ID = (-2)*ones(1,N) ; % init ID
% set INC with art or artm for domain boundary sides (i j)
for r = 1:M % loop over all DOMAIN triangles
tri = TOP(r,:);
tri = [tri tri];
for jj = 1:3 % jj = # of actual triangle side
i = tri(jj+1); j = tri(jj+2); % (i j) = side jj of element r
w = inel(i,j); % find n triangles in TOPD with side (i j) or side (j,i)
if w{1} == 1 % n = w{1} = 1 => (i j) = BOUNDARY side of domain
% (i j) = boundary side
ind = find(BART(:,1) == i);
if isempty(ind) bart = BART(1,2); % default boundary type
else bart = BART(ind(1),2); % extra boundary type
end
if TEST
disp(['boundary side [r i j jj] = ' num2str([r i j jj]) ' is of type = ' num2str(bart)])
end
INC(r,jj) = bart; % jj = boundary side of kind = bart
if bart == 0 % i and j are dirichlet knotes
ID([i j]) = 1;
else % i and j are not dirichlet boundary knotes
if ID(i) ~= 1 ID(i) = 0; end
if ID(j) ~= 1 ID(j) = 0; end
end
end; end; end
if TEST
ind = find(ID == 1)
disp(['all dirichlet knotes']); pause
end