FEM continued, Poisson equ, PDEtool

• [EEHJ]
• [Joh] Claes Johnsson Intro to FEM
• PDE toolbox manual+handouts (Available in math library also)
• A document of PDE Toolbox by Jouni Savolainen (-98)
• [Coo] Cooper

Scanned transparencies

Example PDEtool-session

use matlab
matlab
>> pdetool

• Option: Snap the grid
• Draw unit circle
• Boundary (Dirichlet red, double click on bdary)
• PDE - edit coeffs.
• mesh, refine
• = for solve
• Export everything

>> whos
  Name      Size         Bytes  Class

  a         1x3              6  char array
  ans       1x12            24  char array
  b        10x4            320  double array
  c         1x3              6  char array
  d         1x3              6  char array
  e         7x48          2688  double array
  f         1x3              6  char array
  g        10x4            320  double array
  gd       10x1             80  double array
  ns        3x1             24  double array
  p         2x213         3408  double array
  sf        1x3              6  char array
  t         4x376        12032  double array
  u       213x1           1704  double array

Grand total is 2599 elements using 20630 bytes

>> plot(u)  
>> min(u)

ans =

     0

>> max(u)

ans =

    0.2504

>> pdemesh(p,e,t)
>> pdesurf(p,t,u)
>> help tri2grid

 TRI2GRID Interpolate from triangular mesh to rectangular grid.
 
        UXY=TRI2GRID(P,T,U,X,Y) computes the function values UXY
        over the grid defined by the vectors X and Y, from the
        function U with values on the triangular mesh defined by P and T.
        Values are computed using linear interpolation in
        the triangle containing the grid point. Vectors X and Y
        must be increasing.
        ...
        For grid points outside of the triangular mesh, NaN is
        returned in UXY, TN, A2, and A3.
 
        See INITMESH, REFINEMESH, ASSEMPDE

>> x=linspace(-1,1,20);y=x;
>> uxy=tri2grid(p,t,u,x,y);
>> surf(x,y,uxy)
>> figure
>> surf(x,y,uxy)
>> [X,Y]=meshgrid(x,y);
>> uexa=(1-X.^2-Y.^2)/4;
>> mesh(x,y,uxy-uexa)
>> max(max(uxy-uexa))

ans =

  -6.8863e-06

Basic representation of geometry and grid

• Choose max edge size as 1 (geometry is unit circle).
• Show labels ...
• Initmesh
• Export geometry and grid

------------------------------------
function h=pdeplotn(p,e,t)

h=pdeplot(p,e,t);
hold on


pp = zeros(1,size(p,2)); 
 
for i=1:size(t,2); 
  for j=1:3; 
    tt = t(j,i); 
    if ~pp(tt), 
      text(p(1,tt),p(2,tt),sprintf('%i',tt));
      pp(tt)=1; 
      end
    end
  sum(p(1:2,t(1:3,i))')/3;
  text(ans(1),ans(2),sprintf('%i',i));
  end
hold off;
--------------------

Some more ideas

for i=1:9;Id=eye(9,9);pdesurf(p,t,Id(:,i)); view(30,20);grid;pause;end;

Continue with lecture nr. 8, Apr 13

Shape functions , start with linear

x=[0 1 1];y=[0 0 1];

fill(x,y,'g')   % Triangle with vertices at (x(i),y(i))

z=[1 0 0];fill3(x,y,z,'b');grid;hold on;fill(x,y,'g')

>> view(-30,20)
z=[0 1 0];fill3(x,y,z,'r')
z=[0 0 1];fill3(x,y,z,'y')

% Try various views

>> view(45,0) 
>> view(45,80)
>> view(130,70)

x1=[0 1 0.5];x2=[1 1 0.5];x3=[.5 1 0];x4=[0 .5 0];
y1=[0 0 .5];y2=[0 1 .5];y3=[.5 1 1];y4=[0 .5 1];

>> fill(x1,y1,'y',x2,y2,'b',x3,y3,'g',x4,y4,'r')

z1=[0 0 1];z2=[0 0 1];
fill3(x1,y1,z1,'y',x2,y2,z2,'b');grid
view(30,60) 

Quadratic shape functions:

x=[0  1  2 1.5 1 .5 0];
y=[0 0   0 .5 1 .5 0];
>> fill(x,y,'g')