Программирование МКЭ в МATLAB. Даутов P.З. - 53 стр.

2.3. Вклад элементов в систему МКЭ.                                                                     53


f u n c t i o n [ K, F ] = assemba ( p , t , c , a , b1 , b2 , f )
% ASSEMBA: Assembles a r e a i n t e g r a l c o n t r i b u t i o n s i n a PDE p r o b l e m .
%                  S t i f f n e s s matrix K a s s o c i a t e s with t h e o p e r a t o r
%                            u −−> −d i v ( c grad u ) + b . g r a d u+ au
%
% The f o l l o w i n g c a l l i s a l s o a l l o w e d :
% K = ASSEMBA( p , t , c , a , b1 , b2 )
%
np = s i z e ( p , 2 ) ;
i 1 = t ( 1 , : ) ; i 2 = t ( 2 , : ) ; i 3 = t ( 3 , : ) ; % mesh p o i n t i n d i c e s
s d l= t ( 4 , : ) ; % sub domain l a b e l s

% barycenter of the t r i a n g l e s
x = ( p ( 1 , i 1 )+p ( 1 , i 2 )+p ( 1 , i 3 ) ) / 3 ;
y = ( p ( 2 , i 1 )+p ( 2 , i 2 )+p ( 2 , i 3 ) ) / 3 ;

% g r a d i e n t o f t h e t r i a n g l e ba s e f u n c t i o n s , m u l t i p l i e d on J
g1x=p ( 2 , i 2 )−p ( 2 , i 3 ) ; g1y=p ( 1 , i 3 )−p ( 1 , i 2 ) ; g2x=p ( 2 , i 3 )−p ( 2 , i 1 ) ;
g2y=p ( 1 , i 1 )−p ( 1 , i 3 ) ; g3x=p ( 2 , i 1 )−p ( 2 , i 2 ) ; g3y=p ( 1 , i 2 )−p ( 1 , i 1 ) ;

J=abs ( g 3 y . * g2x−g 3 x . * g2y ) ;         % J=2* a r e a

% e v a l u a t e c , b , a on t r i a n g l e s b a r y c e n t e r
cx     = calc (x , y , sdl , c ) ;
ax     = calc (x , y , sdl , a ) ;
b1x = c a l c ( x , y , s d l , b1 ) ;
b2x = c a l c ( x , y , s d l , b2 ) ;

% d i a g o n a l and o f f d i a g o n a l e l e m e n t s o f mass matrix
ao = ( ax / 2 4 ) . * J ;   ad = 4 * ao ; % 'exa ct ' i n t e g r a t i o n
% ao = ( ax / 1 8 ) . * J ; ad = 3 * ao ; % q u a d r a t u r e r u l e

% c o e f f i c i e n t s o f t h e s t i f f n e s s matrix
cx = ( 0 . 5 * cx ) . /J ;
k12= c x . * ( g 1 x . * g2x+g 1 y . * g2y)+ao ;
k23= c x . * ( g 2 x . * g3x+g 2 y . * g3y)+ao ;
k31= c x . * ( g 3 x . * g1x+g 3 y . * g1y)+ao ;

i f a l l ( b1x==0) && a l l ( b2x==0) % symmetric problem
   K =       s p a r s e ( i 1 , i 2 , k12 , np , np ) ;
   K = K+s p a r s e ( i 2 , i 3 , k23 , np , np ) ;
   K = K+s p a r s e ( i 3 , i 1 , k31 , np , np ) ;

   K = K + K. ' ;

  K = K+s p a r s e ( i 1 , i 1 , ad−k31−k12 , np , np ) ;
  K = K+s p a r s e ( i 2 , i 2 , ad−k12−k23 , np , np ) ;
  K = K+s p a r s e ( i 3 , i 3 , ad−k23−k31 , np , np ) ;
else
  % b contributions
  b1x=b1x / 6 ; b2x=b2x / 6 ;
  bg1=b 1 x . * g1x+b 2 x . * g1y ;
  bg2=b 1 x . * g2x+b 2 x . * g2y ;
  bg3=b 1 x . * g3x+b 2 x . * g3y ;