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

UptoLike

ВУЗ: 

КФУ | Казань

Составители: 

Рубрика: 

2.3. Вклад элементов в систему МКЭ. 51
f u n c t i o n y=c a l c ( x , y , sdl , f )
i f isn u me r ic ( f )
y=f ;
e l s e i f p d eis f u n c ( f ) ,
y=f e v a l ( f , x , y , s d l ) ;
e l s e i f i s c h a r ( f )
y=pdetexpd ( x , y , sd l , [ ] , [ ] , [ ] , [ ] , f ) ;
end
Здесь pdeisfunc(f) и pdetexpd — функции pde toolbox ; pdeisfunc(f)
возвращает 1, если f имя m-файла, mex- или dll-файла в текущем
директории пользователя или в путях поиска функций в Маtlab, а
также, если f встроенная функция MatLab или указатель на функ-
цию. Функция pdetexpd реализует вычисление строк, примеры кото-
рых приведены выше. Отметим, что если входной параметр f есть
число или числовой вектор, то функция calc просто возвращает его.
2.3. Вклад элементов в систему МКЭ.
В разделе 2.1, стр. 45, мы получили все расчетные формулы для
сборки матриц и векторов для P
1
элементов. Сведем в одну функ-
цию вычисление глобальной матрицы жесткости K и вектора сил F,
учитывающих вклад элементов. Начнем с невекторизованной версии
функции.
f u n c t i o n [K, F ] = assemban ( p , t , c , a , b1 , b2 , f )
% ASSEMBAN: Assembles a re a i n t e g r a l c o n t r i b u t i o n s in a PDE p roblem.
% No nv ector iz ed v e r s i o n .
% S t i f f n e s s matrix K a s s o c i a t e s wih the opa ra tor
% u −−> −d iv ( c grad u ) + b. gra d u+ au
%
% The f o l l o w i n g c a l l i s a l s o al lowed :
% K = ASSEMBAN(p , t , c , a , b1 , b2 )
%
np=s i z e (p , 2 ) ; nt=s i z e ( t , 2 ) ;
Kt= z e r o s ( nt , 3 , 3 ) ; % a l l l o c a l s t i f f n e s s matrix
F = z e r o s (np , 1 ) ; % g l o b a l f o r c e v e cto r
f o r i t =1: nt
% mesh poi nt i n d i c e s
i 1=t (1 , i t ) ; i 2=t (2 , i t ) ; i 3=t (3 , i t ) ;
s d l=t (4 , i t ) ; % subdomain l a b e l s
% b a ry c ent e r o f the t r i a n g l e
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 adie n t o f the t r i a n g l e base f unc t i o n s , m u l t i p l i e d on J
2.3. Вклад элементов в систему МКЭ.                                                                 51


f u n c t i o n y=c a l c ( x , y , s d l , f )
i f isnumeric ( f )
    y=f ;
e l s e i f pdeisfunc ( f ) ,
    y=f e v a l ( f , x , y , s d l ) ;
e l s e i f ischar ( f )
    y=pdetexpd ( x , y , s d l , [ ] , [ ] , [ ] , [ ] , f ) ;
end

Здесь pdeisf unc(f ) и pdetexpd — функции pde toolbox ; pdeisf unc(f )
возвращает 1, если f имя m-файла, mex- или dll-файла в текущем
директории пользователя или в путях поиска функций в Маtlab, а
также, если f встроенная функция MatLab или указатель на функ-
цию. Функция pdetexpd реализует вычисление строк, примеры кото-
рых приведены выше. Отметим, что если входной параметр f есть
число или числовой вектор, то функция calc просто возвращает его.

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

   В разделе 2.1, стр. 45, мы получили все расчетные формулы для
сборки матриц и векторов для P1 элементов. Сведем в одну функ-
цию вычисление глобальной матрицы жесткости K и вектора сил F ,
учитывающих вклад элементов. Начнем с невекторизованной версии
функции.
f u n c t i o n [ K, F ] = assemban ( p , t , c , a , b1 , b2 , f )
% ASSEMBAN: 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 .
%                     Nonvectorized v e r s i o n .
%                     S t i f f n e s s matrix K a s s o c i a t e s wih t h e o p a 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 = ASSEMBAN( p , t , c , a , b1 , b2 )
%
np=s i z e ( p , 2 ) ; nt=s i z e ( t , 2 ) ;
Kt= z e r o s ( nt , 3 , 3 ) ; % a l l l o c a l s t i f f n e s s matrix
F = z e r o s ( np , 1 ) ;           % global force vector
f o r i t =1: nt
   % mesh p o i n t i n d i c e s
    i 1=t ( 1 , i t ) ; i 2=t ( 2 , i t ) ; i 3=t ( 3 , i t ) ;
    s d l=t ( 4 , i t ) ; % subdomain l a b e l s

  % barycenter of the t r i a n g l e
  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

Страницы