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

2.4. Учет краевых условий. Формирование системы МКЭ.                                                               59


   s d l=l o c ( eD ( 3 , : ) ) ;

   i 1=eD ( 1 , : ) ; i 2=eD ( 2 , : ) ; % i n d i c e s o f s t a r t and end p o i n t s
   iD =[ i 1 i 2 ] ; % D i r i c h l e t p o i n t s i n d i c e s ( with copy )
   [ id ,∼]= f i n d ( s p a r s e ( iD , iD , 1 , np , np ) ) ; % i d=D i r i c h l e t p o i n t i n d i c e s

   iN=o n e s ( 1 , np ) ; iN ( i d )= z e r o s ( 1 , numel ( i d ) ) ;
   iN=f i n d ( iN ) ;        % i n d i c e s o f non−D i r i c h l e t p o i n t s i n t h e domain
   niN=numel ( iN ) ;
   N=s p a r s e ( iN , 1 : niN , 1 , np , niN ) ;

  i f nargout >=3 % e v a l u a t e uD i n D i r i c h l e t p o i n t s
     uD( i 1 ) = c a l c ( p ( 1 , i 1 ) , p ( 2 , i 1 ) , s d l , bc.uD ) ;
     uD( i 2 ) = c a l c ( p ( 1 , i 2 ) , p ( 2 , i 2 ) , s d l , bc.uD ) ;
  end
end

Здесь мы воспользовались MatLab функцией ismember. Она имеет
несколько вариантов вызова. Команда t = ismember(a, b) возвращает
логический вектор такой же длины как и a, и такой, что tk = 1, если
ak ∈ b (является элементом b) и tk = 0 — иначе.
Формирование и решение системы МКЭ. На основании при-
веденных выше функций можем легко написать различные полез-
ные функции для формирования и решения системы алгебраических
уравнений МКЭ.
    Следующая функция является аналогом функции assempde в pde
toolbox . Она позволяет как находить решение исходной задачи, так и
формировать компоненты системы уравнений МКЭ. Функция легка
для понимания и не требует дополнительных пояснений.
f u n c t i o n [ K, F , N, uD, H,G] = assembpde ( bc , p , e , t , c , a , b1 , b2 , f )
% ASSEMBPDE: Assembles a PDE p r o b l e m .
%
% u=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , f )
% [ K0 , F0 , N, uD]=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , f )
% [ K, F , N, uD , H,G]=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , f )
% [ K0 , M0,N]=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , d )
%
% u=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , f ) a s s e m b l e s and s o l v e s
% t h e PDE problem Lu=−d i v ( c * grad ( u))+ b . g r a d ( u)+a *u=f
% on a mesh d e s c r i b e d by p , e , t , with boundary c o n d i t i o n s
% g i v e n by b c . I t e l i m i n a t e s t h e D i r i c h l e t boundary c o n d i t i o n s
% from t h e system o f l i n e a r e q u a t i o n s when s o l v i n g f o r u .
% The s o l u t i o n u i s r e p r e s e n t e d a s a column v e c t o r o f s o l u t i o n
% v a l u e s a t t h e c o r r e s p o n d i n g node p o i n t s from p .
%
% [ K0 , F0 , N, uD]=ASSEMBPDE( bc , p , e , t , c , a , b1 , b2 , f ) a s s e m b l e s t h e
% PDE problem by e l i m i n a t i n g t h e D i r i c h l e t boundary c o n d i t i o n s