No BSD License  

Highlights from
Two step Taylor Galerkin algorithm

image thumbnail

Two step Taylor Galerkin algorithm

by

 

18 Jun 2007 (Updated )

The two step Taylor Galerkin algorithm applied to nonlinear structural dynamics

Contour1(NN,NE,NEN,X,NOC,FF)
function [Done] = Contour1(NN,NE,NEN,X,NOC,FF)

global hh

% Nodal variable is stored in FF(1:NN)

    FMAX = max(FF); FMIN = min(FF);
    
  NCL = 10; %no. of colours
  STP = (FMAX - FMIN) / NCL;
  % Red-Orange-Green-Blue-Magenta colors 
  % for 10 contour lines [MM(1,:) to MM(10,:)]
  for i=1:10
    a=((10-i)/12); b=1; c=1; 
    H=[a,b,c];MM(i,:)=hsv2rgb(H);
  end;
    
% Find boundary lines
    % Edges defined by nodes in NOC to nodes in NCON
    for IE = 1 : NE
       for I = 1 : NEN
          I1 = I + 1; 
          if I1 > NEN 
             I1 = 1;
          end;
          NCON(IE, I) = NOC(IE, I1);
       end;
    end;
    for IE = 1 : NE
     for I = 1 : NEN
       I1 = NCON(IE, I); I2 = NOC(IE, I);
     	 INDX = 0;
     	 for JE = IE + 1 : NE
     		for J = 1 : NEN
           flow=2;
           if (NCON(JE, J) == 0)
                 flow=1;
           end;
           if ((I1 ~= NCON(JE, J)) & (I1 ~= NOC(JE, J))) 
              flow=1;
           end;
           if ((I2 ~= NCON(JE, J)) & (I2 ~= NOC(JE, J)))
              flow=1;
           end;
           if (flow==2)
              NCON(JE, J) = 0; INDX = INDX + 1;
           end
			end;
       end;
       if INDX > 0
          NCON(IE, I) = 0;
       end;
    end;
   end;
   axis off
   
  %============  Draw Boundary  ==============
  for IE = 1 : NE
     for I = 1 : NEN
        if NCON(IE, I) > 0
           I1 = NCON(IE, I); I2 = NOC(IE, I);
           xe(1)=X(I1,1);xe(2)=X(I2,1);
           ye(1)=X(I1,2);ye(2)=X(I2,2);
           line (xe,ye, 'LineWidth', 2, 'color', 'b')
        end;
     end;
  end
  
  %===========  Contour Plotting  ===========
     for IE = 1 : NE
        if (NEN == 3)
           for IEN = 1 : NEN
              IEE = NOC(IE, IEN);
              U(IEN) = FF(IEE);
              XX(IEN) = X(IEE, 1);
              YY(IEN) = X(IEE, 2);
           end;
           LPLOT(U,XX,YY,FMIN,STP,NCL, MM);
        elseif (NEN == 4)
           XB = 0; YB = 0; UB = 0;
           for IT = 1 : NEN
              NIT = NOC(IE, IT);
              XB = XB + .25 * X(NIT, 1);
              YB = YB + .25 * X(NIT, 2);
              UB = UB + .25 * FF(NIT);
           end;
           for IT = 1 : NEN
              IT1 = IT + 1;
              if (IT1 > 4)
                 IT1 = 1;
              end;              
              XX(1) = XB; YY(1) = YB; U(1) = UB;              
              NIE = NOC(IE, IT);
              XX(2) = X(NIE, 1); YY(2) = X(NIE, 2); U(2) = FF(NIE);
              NIE = NOC(IE, IT1);
              XX(3) = X(NIE, 1); YY(3) = X(NIE, 2); U(3) = FF(NIE);
              LPLOT(U,XX,YY,FMIN,STP,NCL, MM);
           end;
        else
           disp( 'NUMBER OF ELEMENT NODES > 4 IS NOT SUPPORTED')
           stop
        end;
     end;
     
     % Draw legend
        for i=1 : NCL
           cc(i,1)=FMIN+STP*i; 
        end;
        cc1=num2str(cc,3);
        legend(hh,cc1,-1);
        axis equal;
        Done=1;
        
% End of function contour1        
        
%----------------------- function LPLOT             
function [] = LPLOT(U,XX,YY,FMIN,STP,NCL, MM)
global hh
%THREE POINTS IN ASCend;ING ORDER
		     for I = 1 : 2
           C = U(I); II = I;
           for J = I + 1 : 3
              if C > U(J)
                 C = U(J); II = J;
              end;
           end;
           U(II) = U(I); U(I) = C;
           C1 = XX(II); XX(II) = XX(I); XX(I) = C1;
           C1 = YY(II); YY(II) = YY(I); YY(I) = C1;
        end;
        
      SU = (U(1) - FMIN) / STP;
      II = fix(SU+1.e-10)+1;
      if II>NCL
         II=NCL;
      end;
      
       UT = FMIN + II * STP;
        while UT <= U(3)
           X1 = ((U(3)-UT)*XX(1) + (UT-U(1))* XX(3))/(U(3) - U(1));
           Y1 = ((U(3)-UT)*YY(1) + (UT-U(1))* YY(3))/(U(3) - U(1));
           L = 1;
           if UT > U(2) 
             L = 3;
           end;
           X2 = ((U(L)-UT)* XX(2) + (UT-U(2))*XX(L))/(U(L) - U(2));
           Y2 = ((U(L)-UT)* YY(2) + (UT-U(2))*YY(L))/(U(L) - U(2));
           xe(1)=X1;xe(2)=X2;ye(1)=Y1;ye(2)=Y2;           
           hh(II,1)=line(xe,ye,'LineWidth',2,'color',MM(II,:));
           UT = UT + STP;
           II = II + 1;
        end;    
        
% End of function LPLOT
        
       

Contact us