examples.tex

FreeFEM is an implementation of the GFEM language dedicated to the finite element method. It provid

%
% $Id: examples.tex,v 1.4 2001/10/20 23:57:17 prudhomm Exp $
%
% SUMMARY:      
% USAGE:        
%
% AUTHOR:       Christophe Prud'homme <[e-mail]>
% ORG:           MIT
% E-MAIL:       [e-mail]
%
% ORIG-DATE:     8-Feb-97 at 16:47:37
% LAST-MOD: 20-Oct-01 at 15:44:56 by Christophe Prud'homme
%
% DESCRIPTION:  
% This is part of the FreeFEM Documentation Manual
% Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001
%   Christophe Prud'homme and Olivier Pironneau
% See the file fdl.tex for copying conditions.
% DESCRIP-END.


\chapter{Examples}
\label{cha:examples}



%% node  Triangulations examples, Scalar examples, Examples, Examples
\section{Triangulations examples}
%%  node-name,  next,  previous,  up

\begin{itemize}
\item
\textbf{A UNIT RING (INNER RADIUS IS 0.25)}
\begin{Verbatim}
twopi := 2*pi;
border(1,0,twopi,60) { x := cos(t); y := sin(t) };
border(2,0,twopi,20) { x := 0.25*cos(-t); y := 0.25*sin(-t) };
buildmesh(400); 
\end{Verbatim}

\item
\textbf{THE RECTANGLE [(0,0),(0,2),(2,1),(0,10)]}
\begin{Verbatim}
border(1,0,2,20) { x:= t; y:= 0 };
border(1,0,1,10) { x:= 1; y:= t };
border(1,0,2,20) { x:= 2-t; y:= 1 };
border(1,0,1,10) { x:= 0; y:= 1-t };
\end{Verbatim}

\item
\textbf{A SQUARE WITH WELL IDENTIFIED SIDES AND CONTROL OF IB AT CORNERS}
\begin{Verbatim}
border(1,0,4,41)
{
  if(t<=1) then  {  x:=t; y:=0 };
  if((t>1)and(t<2)) then {  x:=1; y:=t-1; ib:= 2 };
  if((t>=2)and(t<=3)) then { x:=3-t; y:=1; ib:= 3 };
  if(t>3) then { x:=0; y:=4-t; ib:= 4 }
};
buildmesh(400);
\end{Verbatim}

\item
\textbf{ MULTI-REGIONS CIRCLE}
\begin{Verbatim}
border(1,0,2,17) {x:= cos(pi*t); y:= sin(pi*t)};
border(0,-1,1,7) { x:= t; y:=0; };
border(0,0,1,4) { x:=0;y:=t };
buildmesh(300);
/* observe the value of "region" by using "show triangle numbers */
\end{Verbatim}
\end{itemize}

%% node  Scalar examples, Complex number example, Triangulations examples, Examples
\section{Scalar examples}
%%  node-name,  next,  previous,  up

\begin{itemize}

\item
\textbf{ ELECTROSTATIC CONDENSOR}
\begin{Verbatim}
/* a circle of radius 5 centered at (0,0) */
border(1,0,2*pi,60) begin x := 5 * cos(t); y := 5 * sin(t) end ;
/* The rectangle on the right */
border(2,0,1,4) begin x:=1+t; y:=3 end ;
border(2,0,1,24) begin x:=2; y:=3-6*t end ;
border(2,0,1,4) begin x:=2-t; y:=-3 end ;
border(2,0,1,24) begin x:=1; y:=-3+6*t end ;
/* The rectangle on the left */
border(3,0,1,4) begin x:=-2+t; y:=3 end ;
border(3,0,1,24) begin x:=-1; y:=3-6*t end ;
border(3,0,1,4) begin x:=-1-t; y:=-3 end ;
border(3,0,1,24) begin x:=-2; y:=-3+6*t end ;
buildmesh(800);

/* Boundary conditions and PDE */
solve(v)
begin
   onbdy(1) v = 0;
   onbdy(2) v = 1;
   onbdy(3) v = -1;
   pde(v) -laplace(v) =0;
end;
plot(v);
\end{Verbatim}

\item
\textbf{HEAT CONDUCTION AND RADIATION }
\begin{Verbatim}
border(1,0,22,89)
begin 
 if(t<=10)then begin x:= t; y:=0 ; ib:=3 end;
 if((t>10)and(t<11))then begin x:=10; y:=t-10; ib:=2 end;
 if((t>=11)and(t<=21))then begin x:=21-t; y:=1; ib:=4  end;
 if(t>21)then begin x:=0; y:=22-t end;
end;
buildmesh(800);

changewait; 
t0 := 10;  t1 := 100;  te := 25; b=0.1; c = 5.0e-8;
w =  (b + 2*c * (te+546)*(te+273)*(te+273));
solve(v,1)
begin
  onbdy(1) v=t0; onbdy(2) v = t1; onbdy(3) dnu(v)=0;
  onbdy(4) id(v) * w + dnu(v) = te * w;
  pde(v) -laplace(v) * y =0;
end;
iter(10)
begin u=v;
  w =  (b + c * (u+te + 546)*((u+273)*(u+273) + (te+273)*(te+273)));
solve(v,-1) begin
  onbdy(1) v=t0; onbdy(2) v = t1;
  onbdy(3) dnu(v)=0; onbdy(4) id(v)*w  + dnu(v)= te * w;
  pde(v) -laplace(v) * y =0; plot(v);
end;
end
\end{Verbatim}

\item
\textbf{HEAT: NON HOMOGENEOUS MATERIAL}
\begin{Verbatim}
r0 := 1.0; r1 := 2.0;
border(1,0,22,89) 
begin  
 region :=1;
 if(t<10)then  begin x:= t; y:=0 ; ib:=3 end;
 if((t>=10)and(t<=11))then  begin x:=10; y:=r1*(t-10); ib:=2 end;
 if((t>11)and(t<21))then begin x:=21-t; y:=r1; ib:=4  end;
 if(t>=21)then begin x:=0; y:=r1*(22-t) end;
end;
border(0,0,10,41) begin x:= t; y:=r0 end;
buildmesh(800);

t0 = 10;  t1 = 100;  te := 25; kappa =0.01 +  max(y-1,0)/(y-1.0001);
solve(v) 
begin
onbdy(1) v=t0; 
onbdy(2) v = t1; 
onbdy(4) dnu(v)=0.2; 
onbdy(3) dnu(v)=0;
pde(v) -laplace(v)*kappa*y +id(v)*kappa*y =0; 
plot(v);
end;
\end{Verbatim}

\item
\textbf{COMPRESSIBLE POTENTIAL FLOW}
\begin{Verbatim}
changewait;/* gamma = 1.4, outer circle radius is 5 */
mach1 := 1/sqrt(6); machinfty = 0.85*mach1; 
rhoinfty=sqrt((1-machinfty^2)^5);
solve(phi) begin
onbdy(1) dnu(phi) = rhoinfty*machinfty*x/5; onbdy(2) dnu(phi) = 0;
pde(phi) id(phi)*0.0001-laplace(phi) = 0;
end;
u1 = dx(phi); u2 = dy(phi); rho=sqrt((1-(u1^2+ u2^2))^5); plot(phi);

iter(5) 
begin 
 solve(phi)
   onbdy(1)  dnu(phi) =rhoinfty*machinfty*x/5; onbdy(2) dnu(phi) = 0;
   pde(phi)  id(phi)*0.0001-laplace(phi)*rhop=0;
 end;
 u1 = dx(phi); u2 = dy(phi); rho=sqrt((1-(u1^2+ u2^2))^5);
 rhop = convect(rho,u1,u2,0.1); plot(rho)
end;
plot(sqrt((u1^2+u2^2))/mach1);
\end{Verbatim}

\item
\textbf{NAVIER STOKES EQUATIONS}
\begin{Verbatim}
/* Poor but better than none algorithm*/
changewait;
border(1,0,1,6)  begin x:=0;      y:=1-t   end;
border(2,0,1,15) begin x:=2*t;    y:=0     end;
border(2,0,1,10) begin x:=2;      y:=-t    end;
border(2,0,1,20) begin x:=2+3*t;  y:=-1    end;
border(2,0,1,35) begin x:=5+15*t; y:=-1    end;
border(3,0,1,10) begin x:=20;     y:=-1+2*t end;
border(4,0,1,35) begin x:=5+15*(1-t); y:=1 end;
border(4,0,1,40) begin x:=5*(1-t);y:=1     end;
buildmesh(900);

nu = 0.002; dt := 0.4;

/* initial pressure */
        solve(p,1) 
        onbdy(1)dnu(p) =-2*nu; 
        onbdy(3) p=0; onbdy(2,4) dnu(p) = 0;
        pde(p) - laplace(p)= 0;
        end;
/* initial horizontal velocity */
        solve(u,2) begin
        onbdy(1) u = y*(1-y);
        onbdy(3) dnu(u) = 0; onbdy(2,4) u = 0;
    pde(u) id(u)/dt-laplace(u)*nu = -min(y*y-y,0)/dt;
    end;        
/* initial vertical velocity */
        solve(v,3)begin
        onbdy(1,3)v = 0; onbdy(2,4) v = 0;
    pde(v) id(v)/dt-laplace(v)*nu =0; 
    end;
    un = u; vn = v; 
iter(80)
begin f=convect(un,u,v,dt);  g=convect(vn,u,v,dt);
/*Horizontal velocity*/
    solve(u,-2) begin
    onbdy(1) u = y*(1-y); onbdy(2,4) u = 0; 
    onbdy(3)dnu(u)=0; 
    pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
    end;
    plot(u);
/* Vertical velocity */
        solve(v,-3) begin
        onbdy(1,2,3,4) v = 0;
    pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
    end;
/*  Pressure */
        solve(p,-1) begin
        onbdy(1)dnu(p) =-2*nu;
        onbdy(3) p=0; onbdy(2,4) dnu(p) = 0;
        pde(p) -laplace(p)= -(dx(f) + dy(g))/dt;
        end;
    un = u; vn = v;
 end ; 
    save('u.dta',u); save('v.dta',v); save('p.dta',p); plot(u);
\end{Verbatim}
\end{itemize}

%% node  Complex number example, 2-system example, Scalar examples, Examples
\section{Complex number example}
%%  node-name,  next,  previous,  up

\begin{Verbatim}
complex; nowait;
border(1,0,1,10) begin x:=t; y:=0; end;
border(1,0,1,10) begin x:=1; y:=t; end;
border(2,0,1,10) begin x:=1-t; y:=1; end;
border(1,0,1,10) begin x:=0; y:=1-t; end;
buildmesh(200);

solve(u)  /* observe than Re(u) = Im(u) */
begin
     onbdy(1,2) u=0;
      pde(u) id(u)-laplace(u)=1+I ;
end;
v=Im(u);
 plot(u);plot(v);plot(u-v);
\end{Verbatim}

%% node  2-system example,  , Complex number example, Examples
\section{2-system example}
%%  node-name,  next,  previous,  up

\begin{Verbatim}
/* This is a 2-system example for which the solution is know
analytically, thus the precision of Gfem can be checked */
nowait;
ns:=40; 
border(1,0,2*pi,2*ns) begin x:= 3*cos(t); y:= 2*sin(t); end;
border(2,0,2*pi,ns)   begin x:= cos(-t); y:= sin(-t); end;
buildmesh(ns*ns);

ue= sin(x+y);
ve = ue;
p = ue;
nx = -x;
ny =- y;
dxue = cos(x+y);

c = 0.2;
a1 = y;
a2 =x;
nu = 1;
nu11 = 1; 
nu22 = 2;
nu21 =0.3;
nu12 =0.4;
b=1;

dnuue=dxue*(nu*(nx+ny) + 
(nu11 + nu12)*nx + (nu21+ nu22)*ny);
g = ue*c+dnuue;
f = b*ue+dxue*(a1+a2) +ue*(2*nu+nu11+nu12+nu21+nu22);

solve(u,v) begin
onbdy(1) u = p;
onbdy(1) v = p;
onbdy(2) 
          id(u)*c/2 + id(v)*c/2 + dnu(u) = g;
onbdy(2)   id(v)*c + dnu(v) = g;
pde(u) id(u)*b + dx(u)*a1 + dy(u)*a2
          -laplace(u)*nu - dxx(u)*nu11 -   
             dxy(u)*nu12 - dyx(u)*nu21 - dyy(u)*nu22 =f;

pde(v) id(v)*b/2+id(u)*b/2 
       + dx(v)*a1 + dy(v)*a2 -laplace(v)*nu 
        - dxx(v)*nu11 - dxy(v)*nu12 -          
        dyx(v)*nu21 - dyy(v)*nu22 =f;
end;

plot(abs(u-ue) + abs(v-ve));
\end{Verbatim}
%
%%%%%%%%%%%%% Some Settings for emacs and auc-TeX
% Local Variables:
% TeX-master:  "freefem"
% TeX-parse-self: t
% TeX-auto-save: t
% TeX-auto-regexp-list: TeX-auto-full-regexp-list
% End:
%