Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
1d heat finite element

Subject: 1d heat finite element

From: RICARDO

Date: 10 Dec, 2012 21:49:22

Message: 1 of 2

Hello forum.
I am solving the heat equation i 1d using FEM. It works fine for initial condition. However when I tried to run it for t>0 it does not work properly.

Here is my code, you can run it with this inputs and see
[Uh]=ParaFEM1D(1,0,1,10,7)


function [Uh, xnod]=ParaFEM1D(Tf,a, b, Nt, xnel)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This function shows how to set up a finite element solution
%% of a two-point BVP
%% -(ku')' =f, u(a)=ua, u(b)=ub
%% using Galerkin linear Finite Elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function sol = solfem (a,b,xnel)
%% k,f,ua,ub are provided as functions in code
%% numerical integration is used in computation of stiffness matrix and rhs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT:
%% <a,b> = endpoints of interval
%% <xnel> = the number or location of elements
%% if xnel is a scalar, then xnel = number of elements,
%% if xnel is a vector, then it contains position of first node of each element
%% OUTPUT:
%% <sol> is the computed numerical solution at nodes
%% <xnod> position of all nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
set(0,'DefaultLineLineWidth',3);
%%%%%%%%%%%%%%%%%%%%%%%%% set up grid
dt = Tf/Nt;
t = 0:dt:Tf;

if size(xnel) == 1 %% uniform grid
    
   nels = xnel;
   hel = (b-a)./nels;
   xel = a : hel : b-hel;

elseif size(xnel,1) == 1

   nels = size(xnel,2);
   xel=xnel;
   
   for i=1:nels-1
       
       hel(i) = xnel(i+1)-xnel(i);
   
   end;

else

   error('wrong number of elements/grid');

end;

%% set up uniform order of elements = 1 + degree of polynomial =
%% = number of degrees of freedom

ord = zeros(nels,1) + 2; %% type of elements: linear
maxord = max(ord);

%% number of nodes
nnodes = sum(ord-1)+1;

%% derive global indexing of nodes:
%% nod(i,j) is the global number of j'th node in element i
nod = zeros(nels,maxord); myel = zeros(nnodes,2);
n = 1;
for i = 1:nels
    for j = 1:ord(i)
        nod(i,j) = n;
        if j == 1
            myel(n,2) = i;
        elseif j == ord(i)
            myel(n,1) = i;
        else
            myel(n,1) = i;
            myel(n,2) = i;
        end;
        
        if j ~= ord(i)
            n = n+1;
        end
    end;
end;

%% xnod (i=1..nnodes): coordinates of node i
xnod = zeros(nnodes,1);
for i=1:nels-1
    h = xel(i+1)-xel(i);
    hi = h/(ord(i)-1);
    for j=1:ord(i)
        xnod (nod(i,j)) = xel(i) + hi*(j-1);
    end;
end;

i = nels;
h = b-xel(i);
hi=h/(ord(i)-1);

for j=1:ord(i)
    xnod (nod(i,j)) = xel(i) + hi*(j-1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%% set up numerical integration

%% set up quadrature parameters on the reference element (-1,1)
%% set up number of integration points nw, nodes xw, and weights w

if maxord == 1 %% exact for linears
      nw = 1;
      xw(1) = 0.;
      w(1) = 2.;
elseif maxord == 2 %% exact for cubics
      nw = 2;
      xw(1) = -1/sqrt(3); xw(2) = -xw(1);
      w(1) = 1; w(2) = 1;
elseif maxord == 3 %% exact for polynomials of degree 5
      nw = 3;
      xw(1) = -sqrt(3./5.); xw(2)=0.; xw(3) =- xw(1);
      w(1) = 5./9.; w(2)=8./9.; w(3)=w(1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%% matrix and rhs of linear system
b_global = sparse(nnodes,nnodes); m_global = sparse(nnodes,nnodes);
gs_global = sparse(nnodes,nnodes); rhsf = zeros(nnodes,1);
for el = 1:nels
    x1 = xnod(nod(el,1)); %% left endpoint
    x2 = xnod(nod(el,2)); %% right endpoint
    dx = (x2-x1)/2.; %% Jacobian of transformation
    %% compute element stiffness matrix and load vector
    b_local = zeros(ord(el),ord(el)); %% element stiffness matrix
% gs_local = zeros(ord(el),ord(el));
    m_local = zeros(ord(el),ord(el));
% gs_local = zeros(ord(el),ord(el)); %% GS = B+dt*A
    f = zeros(ord(el),1); %% element load vector
    %nw is the order of Gaussian Integration
    for i = 1:nw
        x = x1 + (1 + xw(i))*dx; %% x runs in true element,
                                                %% xw runs in reference element
        [psi,dpsi] = shape(xw(i),ord); %% calculations on ref.element
% kval = feval(@kfun,x);
% fval = feval(@rhsfun,0,x);
        fval = feval(@exfun,0,x);
        f = f + fval * psi * w(i)*dx;
        
        b_local = b_local + (psi*psi')* w(i)*dx; %mass vector
        m_local = m_local + w(i)*(dt*(dpsi*dpsi')/dx/dx) * dx;%stiffness matrix
% gs_local = gs_local + w(i)*(dt*(dpsi*dpsi')/dx/dx) * dx + (psi*psi')* w(i)*dx ;
     end
    %% add the computed element stiffness matrix and load vector to
    %% the global matrix and vector
    rhsf(nod(el,:)) = rhsf(nod(el,:)) + f(:); % F = load vector
    b_global(nod(el,:),nod(el,:)) = b_global(nod(el,:),nod(el,:)) + b_local; % B = mass matrix
    m_global(nod(el,:),nod(el,:)) = m_global(nod(el,:),nod(el,:)) + m_local; % A = stiffness matrix
% gs_global(nod(el,:),nod(el,:)) = gs_global(nod(el,:),nod(el,:)) + gs_local;
end

%% impose Dirichlet boundary conditions
xa = exfun(0,a); xb = exfun(0,b) ;
Uh = zeros(Nt,nnodes);
for i=1:nnodes
% rhsf(i) = rhsf(i) - m_global(i,1)*xa - m_global(i,nnodes)*xb ;
     rhsf(i) = rhsf(i) - (dt*m_global(i,1)-b_global(i,1))*xa - ...
         (dt*m_global(i,nnodes)-b_global(i,nnodes))*xb ;
    m_global(i,1) =0; m_global(1,i)=0;
    b_global(i,1) =0; b_global(1,i)=0;
    m_global(i,nnodes)=0; m_global(nnodes,i)=0;
    b_global(i,nnodes)=0; b_global(nnodes,i)=0;
end

m_global(1,1) = 1; b_global(1,1) = 1; rhsf(1,1) = xa;
m_global(nnodes,nnodes)=1;b_global(nnodes,nnodes)=1; rhsf(nnodes,1) = xb;

for n=1:Nt
   fprintf('itetation %d of %d\n',n,Nt);
% k = t(n+1) - t(n); % time step;
   if n == 1
      Uh(n,:) = b_global\rhsf;
   else
% Uh(ntime,:) = Uh(ntime-1,:)*(gs_global\(dt*b_global));
      Uh(n,:) = (b_global+dt*m_global)\(b_global*Uh(n,:)'+dt*rhsf);
   end
  figure(1)
  exfun(dt*(n-1),xnod)'
  plot(xnod,exfun(dt*(n-1),xnod),'k',xnod,Uh(n,:),'r--');
  pause(0.5),legend('Exact solution','Numerical solution (linear FE)');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% end of algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y,dy] = shape(x,n)

%% shape function on reference element (-1,1)
%% n = 2: linear
%% n = 3: quadratic (must be coded)

if n == 2

   y (1,:) = .5.*(1-x);
   y (2,:) = .5.*(1+x);
   dy (1,:) = -.5;
   dy (2,:) = .5;

end


function y = exfun(t, x)
     
y= sin(pi*x).*exp(-pi^2*t) + sin(2*pi*x).*exp(-4*pi^2*t);
% dy = pi*cos(pi*x).*exp(-pi^2*t) + 2*pi*cos(2*pi*x).*exp(-4*pi^2*t);

Please help

Subject: 1d heat finite element

From: Nasser M. Abbasi

Date: 10 Dec, 2012 22:13:35

Message: 2 of 2

On 12/10/2012 3:49 PM, RICARDO wrote:
> Hello forum.
> I am solving the heat equation i 1d using FEM. It works fine for
>initial condition. However when I tried to run it for t>0 it does not work properly.
>


how about first fixing these issues which Matlab tells you about?

25: The value assigned to variable 't' might be unused.
31: Variable 'hel' might be set by a nonscalar operator.
31: Variable 'hel' might be set by a nonscalar operator.
40: The value assigned to variable 'hel' might be unused.
40: The variable 'hel' appears to change size on every loop
iteration. Consider preallocating for speed.
120: The value assigned to variable 'gs_global' might be unused.

and few more. You can see these from clicking on tools->code analyzer.

I am not saying these will fix or not fix, your logic errors, but
it won't hurt to fix these first.

I think you are asking here for someone to debug your FEM program?
I hope someone will have to do that. FEM is interesting subject.

If I were you, I'd make a print out of the code, and sit
and go over it carefully to make sure it is doing what
it is supposed to do. debugging is an important skill that needs
to be learned in order to become a good programmer.

--Nasser

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us