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

§ 5. P2 сетки                                                                                                         35


np=s i z e ( p , 2 ) ; ne=s i z e ( e , 2 ) ; nt=s i z e ( t , 2 ) ;

%−−−−− s e t e l e m e n t s o f edge o r i e n t e d mesh data −−−−−−−−−−−−
% m a t r i c e s ep , te , e t a r e t h e same a s i n f u n c t i o n adjmeshP1
%−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
% s e t ep = f i r s t 2 rows o f e e
t t =( t ( 1 : 3 , : ) ) ' ;
ep = [ t t ( : , [ 1 , 2 ] ) ; t t ( : , [ 2 , 3 ] ) ; t t ( : , [ 3 , 1 ] ) ] ; % a l l e d g e s − not unique
 [ ep ,∼ , j ]= unique ( s o r t ( ep , 2 ) , 'rows' ) ;                         % unique mesh e d g e s
ep=ep ' ; j=j ' ; nep=s i z e ( ep , 2 ) ; % ep i s s o r t e d : 1 row and each column

% set te
mt = 1 : nt ;
t e = [ j (mt ) ; j (mt+nt ) ; j (mt+2* nt ) ] ; % e l e m e n t s e d g e s

% s e t midpoint i n d i c e s i n e l e m e n t s c o n n e c t i v i t y matrix t 2 .
t 2 =[ t ( 1 : 3 , : ) ; np+t e ; t ( 4 , : ) ] ; % midpoint <−−> edge

% s e t e t = 7 , 8 rows o f e e ( f i r s t attempt )
e t =z e r o s ( 2 , nep ) ; % e d g e 2 e l e m e n t c o n n e c t i v i t y matrix
f o r k=1: nt
    f o r j =1:3
        i t=t e ( j , k ) ;
        if      e t ( 1 , i t )==0, e t ( 1 , i t )=k ;
        e l s e e t ( 2 , i t )=k ; end
    end
end

% s o r t f i r s t 2 rows o f e a s i n ep
f o r j =1: ne
    i f e ( 1 , j )>e ( 2 , j ) ;
       c=e ( 1 , j ) ; e ( 1 , j )=e ( 2 , j ) ; e ( 2 , j )=c ;
       c=e ( 3 , j ) ; e ( 3 , j )=e ( 4 , j ) ; e ( 4 , j )=c ;
       c=e ( 6 , j ) ; e ( 6 , j )=e ( 7 , j ) ; e ( 7 , j )=c ;
    end
end
[∼ , ib , i t ] = i n t e r s e c t ( ep ' , e ( 1 : 2 , : ) ' , 'rows' ) ;

% s e t midpoint i n d i c e s i n boundary c o n n e c t i v i t y matrix e 2 .
e2 =[ e ( 1 : 2 , : ) ; z e r o s ( 1 , ne ) ; e ( 3 : 7 , : ) ] ;
te = [ te ; t ( 4 , : ) ] ;

% s e t 3 : 6 rows o f e e
e e =[ ep ; np +(1: nep ) ; z e r o s ( 5 , nep ) ] ; e e ( 4 : 6 , i b )=e2 ( 4 : 6 , i t ) ;

e2 ( 3 , i t )= e e ( 3 , i b ) ;

% s e t 7 , 8 rows o f e e
f o r i i =1: nep % i i =edge
    i t 1=e t ( 1 , i i ) ; i t 2=e t ( 2 , i i ) ; % n e i g h b o u r s t r i a n g l e s
    i=ep ( 1 , i i ) ; j=ep ( 2 , i i ) ; k= s e t d i f f ( t t ( i t 1 , : ) , [ i j ] ) ; % t 1 v o r t e x i n d i c e s
   % oriented area of t r i a n g l e ( i , j , k)
    d12 = p ( : , j )−p ( : , i ) ;
    d13 = p ( : , k)−p ( : , i ) ;