% 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