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

52                                                Глава 3. Программирование сборки матриц МКЭ


     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 ; % 'exact ' 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

   % 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 ;

   % 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 ;

  Kt ( i t , : , : ) = [ ad−k31−k12+bg1               k12+bg2                        k31+bg3
                         k12+bg1                      ad−k12−k23+bg2                 k23+bg3
                         k31+bg1                      k23+bg2                        ad−k23−k31+bg3 ] ;
  i f n a r g o u t==2
     fx = calc (x , y , sdl , f ) ;
     I =[ i 1 ; i 2 ; i 3 ] ;
     F( I ) = F( I )+ f x * J / 6 ;
  end
end

K=s p a r s e ( np , np ) ; % g l o b a l s t i f f n e s s matrix
f o r i =1:3 , f o r j =1:3
      K=K+s p a r s e ( t ( i , : ) , t ( j , : ) , Kt ( : , i , j ) , np , np ) ;
end , end

i f n a r g o u t ==1, F = [ ] ; end

   Простой анализ этой функции показывает, что для векторизации
достаточно убрать цикл по it, векторизовать команды помеченные
знаком процента с цифрами 1, 2, а также рассылать элементы сразу,
как только их вычислили. В результате получим следующую векто-
ризованную версию функции.