Chapter 4  Examples

4.1  Triangulations examples

• A UNIT RING (INNER RADIUS IS 0.25)

  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);

• THE RECTANGLE [(0,0),(0,2),(2,1),(0,10)]

  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 };

• A SQUARE WITH WELL IDENTIFIED SIDES AND CONTROL OF IB AT CORNERS

  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);

• MULTI-REGIONS CIRCLE

  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 */

4.2  Scalar examples

• ELECTROSTATIC CONDENSOR

  /* 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);

• HEAT CONDUCTION AND RADIATION

  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

• HEAT: NON HOMOGENEOUS MATERIAL

  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;

• COMPRESSIBLE POTENTIAL FLOW

  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);

• NAVIER STOKES EQUATIONS

  /* 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);

4.3  Complex number example

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);

4.4  2-system example

/* 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));