Matlab - Programa de metodo elemntos finitos libro Young W Kwon

% Exemplo 351.m

% para resolver a EDO dad por:

% au''+bu'+cu=1, 0<x<1

% u(0)= e u(1)=0

% usando 5 elementos

%

% descripcao de variaveis

% k= matriz elemento

% f= vetor elemento

% kk= matriz sistema

% ff= vetor sistema

% index= um vetor contendo um sistema associado dofs de cada elemento

% bcdof= um vetor contendo dofs associado com condicoes de fronteira

% bcval= um vetor contendo condicao de valores na fronteira associado com o dofs em bcdof

%

%--------------------------------------------------------------------------%

% entrada de datos

%-------------------------------------------------------------------------

nel=5;    % numero de elementos

nnel=2;    % numero de nodos por elemento

ndof=1;    % numero de dofs por nodo

nnode=6;     %numero total de nodos del sistema

sdof=nnode*ndof;    %dofs do sistema total

%

%--------------------------------------------

% datos de entrada de valores coordenados de cada nodo

%-------------------------------------------

gcoord(1)=0.0 ; gcoord(2)=0.2; gcoord(3)=0.4; gcoord(4)=0.6;

gcoord(5)=0.8; gcoord(6)=1.0;

%

%-------------------------------------------

% datos de entrada de conectividade nodal para cada elemento

%-------------------------------------------

nodes(1,1)=1; nodes(1,2)=2; nodes(2,1)=2; nodes(2,2)=3;

nodes(3,1)=3; nodes(3,2)=4; nodes(4,1)=4; nodes(4,2)=5;

nodes(5,1)=5; nodes(5,2)=6;

%

%------------------------------------------------------

% coeficientes de entrada da EDO

%-----------------------------------------------------

acoef=1;

bcoef=-3;

ccoef= 2;

%

%-----------------------------------------

%datos de entrada para a condicao de fronteira

%-------------------------------------------

bcdof(1)=1;

bcval(1)=0;

bcdof(2)=6;

bcval(2)=0;

%

%------------------------------------------------

% inicializacao de vetores e matrizes

%-----------------------------------------------

ff= zeros(sdof,1);     % inicializacao do sistema do vetor forca

kk=zeros(sdof,sdof);           % inicializacao do sistema matricial

index=zeros(nnel*ndof,1);           % inicializacao do vetor index

%

%------------------------------------------------

% calculo de elemntos das matrizes e vetores e seu emsamble

%---------------------------------------------------

for iel=1:nel                       % loop para o numero total de elementos

%

nl=nodes(iel,1); nr=nodes(iel,2);

xl=gcoord(nl); xr=gcoord(nr);

eleng=xr-xl;

index=feeldof1(iel,nnel,ndof);

%

k= feode2l(acoef,bcoef,ccoef,eleng);

f= fef1l(xl,xr);

[kk,ff]=feasmbl2(kk,ff,k,f,index);

%

end

%

%----------------------------------------

%aplicando condicoes de fronteira

%-----------------------------------------

[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);

%kk, ff, k, f, index

%--------------------------------------

% resolve a equacao matricial

%------------------------------------

fsol=kk\ff;

%

%------------------------------------------------

c1= 0.5/exp(1);

c2=-0.5*(1+1/exp(1));

for i=1:nnode

    x=gcoord(i);

    esol(i)=c1*exp(2*x)+c2*exp(x)+1/2;

end

%

%------------------------------------------

% imprime ambas solucoes exata e por mef

%----------------------------------------

num=1:1:sdof;

results=[num' fsol  esol']

%-----------------------------------------------

function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)

%--------------------------------------------------------

% aqui ha dof=2 e 10

%bcdof(1)=2 e bcdof(2)=10;

%bcval(1)=1.0 e bcval(2)=2.5

%---------------------------------------------

n=length(bcdof);

sdof=size(kk);

%

for i=1:n

    c=bcdof(i);

    for j=1:sdof

        kk(c,j)=0;

    end

    %

    kk(c,c)=1;

    ff(c)=bcval(i);

end

end

function [index]=feeldof1(iel,nnel,ndof)

edof=nnel*ndof;

start=(iel-1)*(nnel-1)*ndof;

%

for i=1:edof

    index(i)=start+i;

end

end

function [f]=fef1l(xl,xr)

%----------------------------------

% elemnto vetor para f(x)=1

% usando elemento linear

 eleng=xr-xl;

f=[eleng/2;eleng/2];

function [k]=feode2l(acoef,bcoef,ccoef,eleng)

%---------------------------------------

a1=-(acoef/eleng); a2=bcoef/2; a3=ccoef*eleng/6;

k=[a1-a2+2*a3,  -a1+a2+a3;...

    -a1-a2+a3 , a1+a2+2*a3];

%----------------------------

function [kk,ff]=feasmbl2(kk,ff,k,f,index)

%--------------------------

edof=length(index);

for i=1:edof

    ii=index(i);

    ff(ii)=f(ii)+f(i);

    for j=1:edof

        jj=index(j);

        kk(ii,jj)=kk(ii,jj)+k(i,j);

    end

end

function [kk,ff]=feasmbl2(kk,ff,k,f,index)

%--------------------------

edof=length(index);

for i=1:edof

    ii=index(i);

    ff(ii)=ff(ii)+f(i)   ;

    for j=1:edof

        jj=index(j);

        kk(ii,jj)=kk(ii,jj)+k(i,j);

    end

    end

end

% -----------------------------------------------

>> metodos_elementos_finitos_young

results =

    1.0000   -0.0000         0

    2.0000   -0.0621   -0.0610

    3.0000   -0.1133   -0.1110

    4.0000   -0.1388   -0.1355

    5.0000   -0.1142   -0.1111

    6.0000         0   -0.0000