Gfem does three tasks:
Gfem is driven by programs read from disk or written in Edit windows.
In integrated environments like MacGfem, PCGfem or XGfem, to execute a program you may use the menu Run or select part of the text and press ENTER(not return). Without menu-window environment, the name of the program file is read at launch time.
This simple program:
border(1,0,6.28,20) begin x:=cos(t); y:=sin(t) end; buildmesh(300); f=x*y; plot(f);
will create a triangulation for the unit circle with 20 points on the boundary , which has received label 1, and a maximum of 300 vertices, then display the level lines of function f=x*y;
If to the previous program
solve(u) begin onbdy(1) u=0; pde(u) -laplace(u)=f end;
is added then the solution is generated for
-laplace(u)=f in Omega , u=0 on the boundary of Omega
Gfem has a small language which generally follows the syntax of the language Pascal. See the list below for the reserved words of this language.
The reserved word begin can be replaced by { and end }
C programmers: caution the syntax is "...};" while most C constructs use " ...;}"
Example 1: Triangulates a circle and plot f = x*y:
border(1,0,6.28,20) begin x:=cos(t); y:=sin(t); end; buildmesh(200); f=x*y; plot(f);
Example 2: on the circle solve the Dirichlet problem -laplace(u)= x*y with u=0 on the boundary:
border(1,0,6.28,20) begin x:=cos(t); y:=sin(t); end; buildmesh(200); solve(u) begin onbdy(1) u =0; pde(u) -laplace(u) = x*y ; end; plot(u);
begin, end, {, },
x, y, t, ib, iv, region, nx, ny,
ln, exp, sqrt, abs, one, min, max,
sin, asin, tan, atan, cos, acos,
cosh, sinh, tanh, pi,
I, Re, Im, complex, scal
dx, dy, dxx, dyy, dxy, dyx, convect,
if, then, else, or, and,iter,
bdy, buildmesh, savemesh, loadmesh,
solve, pde, id, dnu, laplace, div,
onbdy, plot, plot3d, save, load,
halt, include, evalfct, exec,saveall,
user, precise, wait, nowait, changewait,
adaptmesh.
To create a triangulation you must either
In integrated environments, once created, triangulations can be displayed, stored on disk in Gfem format or in text format or even a zoom of its graphic display can be stored in Postscript format (to include it in a TeX file for example).
Gfem stores for each vertex its number and boundary identification number and for each triangle its number and region identification number. Edges number are not stored, but can be recovered by program if needed.
Use it to triangulate domain defined by its boundary. The syntax is
border(ib,t_min,t_max,nb_t) begin ...x:=f(t); ...y:=g(t)... end; buildmesh(nb_max);
where each line with border could be replaced by a line with polygon
polygon(ib,'file_name'[,nb_t]);
where f,g are generic functions and the [...] denotes an optional addition. The boundary is given in parametric form. The name of the parameter must be t and the two coordinates must be x and y. When the parameter goes from t_min to t_max the boundary must be scanned so as to have Omega on its left, meaning counter clockwise if it is the outer boundary and clockwise if it is the boundary of a hole. Boundaries must be closed but they may be entered by parts, each part with one instruction border, and have inner dividing curves; nb_t points are put on the boundary with values t = t_min + i * (t_max - t_min) / (nb_t-1) where i takes integer values from 0 to nb_t-1. The triangulation is created by a Delaunay-Voronoi method with nb_max vertices at most. The size of the triangles is monitored by the size of the nearest boundary edge. Therefore local refinement can be achieved by adding inner artificial boundaries.
Triangulation may have boundaries with several connected components. Each connected component is numbered by the integer ib.
Inner boundaries (i.e. boundaries having the domain on both sides) can be useful either to separate regions or to monitor the triangulation by forcing vertices and edges in it. They must be oriented so that they leave Omega on their right if they are closed. If they do not carry any boundary conditions they should be given identification number ib=0. The usual if... then ... else statement can be used with the compound statement: begin...end. This allows piecewise parametric definitions of complicated or polygonal boundaries.
The boundary identification number ib can be overwritten.For example:
border(2,0,4,41) begin
if(t<=1)then { x:=t; y:=0 };
if((t>1)and(t<2))then { x:=1; y:=t-1; ib=1 };
if((t>=2)and(t<=3))then { x:=3-t; y:=1 };
if(t>3)then { x:=0; y:=4-t }
end;
buildmesh(400);
Recall that begin and { is the same and so is end and }. Here one side of the unit square has ib=1. The 3 other sides have ib=2.
The keyword polygon causes a sequence of boundary points to be read from the file file_name which must be in the same drectory as the program. All points must be in sequential order and describing part of the boundary counter clockwise; the last one should be equal to the first one for a closed curve.
The format is
x[0] y[0] x[1] y[1] ....
each being separated by a tab or a carriage return. The last parameter nb_t is optional; it means that each segment will be divided into nb_t+1 equal segments (i.e. nb_t points are added on each segments).
For example
polygon(1,'mypoints.dta',2); buildmesh(100);
with the file mypoints.dta containing
0. 0. 1. 0. 1. 1. 0. 1. 0. 0.
triangulates the unit square with 4 points on each side and gives ib=1 to its boundary. Note that polygon(1,'mypoints.dta') is like polygon(1,'mypoints.dta',0).
Inner boundaries which divide the domain into several simply connected components are useful to define piecewise discontinuous functions such as dielectric constants, Young modulus...
Inner boundaries may meet other boundaries only at their vertices. Such inner boundaries will split the domain in several sub-domains.
A sub-domain is a piece of the domain which is simply connected and delimited by boundaries.
Each sub-domain has a region number assigned to it. This is done by Gfem, not by the user. Every time border is called, an internal number ngt is incremented by 1. Then when the key word border is invoked the last edge of this portion of boundary assigns this number to the triangle which lies on its left. Finally all triangles which are in this subdomain are reassigned this number.
At the end of the triangulation process, each triangle has a well defined number ngt. The number region is a piecewise linear continuous interpolation at the vertices of ngt. To be exact, the value of region at a vertex (x_0,y_0) is the value of ngt at (x_0,y_0-10^-6) , except if precise is set in which case region is equal to ngt.
Functions are either read or created.
Most usual functions can be used:
max, min, abs, atan, sqrt, cos, sin, tan, acos, asin, one, cosh, sinh, tanh, ln, exp
one(x+y<0) for instance means 1 if x+y<0 and 0 otherwise.
Operators:
and, or, < , <=, < , >=, ==, +, -, *, /, ^ x^2 means x*x
Functions created by a program are displayed only if the key word plot() or plot3d() is used ( here plot(f) ). Derivatives of functions can be created by the keywords dx() and dy().
Unless precise is set, they are interpolated so the results is also continuous piecewise linear (or discontinuous when precise is set). Similarly the convection operator convect(f,u1,u2,dt) defines a new function which is approximately
f(x-u1(x,y)dt , y-u2(x,y)dt) Scalars are also helpful to create functions. Since no data array is attached to a scalar the symbol := is useful to create them, as in
a:= (1+sqrt(5))/2; f= x*cos(a*pi*y);
Here f is a function, a is a scalar and pi is a (predefined) a scalar.
It is possible to evaluate a function at a point as in a:=f(1,0)(see section Value of a function at one point). Here the value of a will be 1 because f(1,0) means f at x=1 and y=0.
There are 6 predefined functions: x,y,ib,region, nx, ny.
The usual if... then ... else statement can be used with an important restriction on the logical expression which must return a scalar value:
if( logical expression) then
{
statement;
....;
statement;
}
else
{
.....
};
The logical expression controls the if by its return being 0 or >0, it is evaluated only once (i.e. with x,y being the coordinates of the first vertex, if there are functions inside the logical expression). Auxiliary variables can be used.
In order to minimize the memory the symbol := tells the compiler not to allocate a data array to this variable. Thus v=sin(a*pi*x); generates an array for v but no array is assigned to a in the statement a:=2.
If f has been defined earlier then it is possible to write a:=f(1.2,3.1); Then a has the value of f at x=1.2 and y=3.1.
It is also allowed to do
x1:=1.2;
y1:=1.6;
a:=f(x1,2*y1);
save("a.dta",a);
Remark: Recall that ,a being a scalar, its value is appended to the file a.dta.
dx(f) is the partial derivative of f with respect to x; the result is piecewise constant when precise is set and interpolated with mass lumping as a the piecewise linear function when precise is not set.
Note that dx() is a non local operator so statements like f=dx(f) would give the wrong answer because the new value for f is place before the end of the use of the old one. The Finite Element Method does not handle convection terms properly when they dominate the viscous terms: upwinding is necessary; convect provides a tool for lagrangian upwinding. By g=convect(f,u,v,dt) Gfem construct an approximation of
f(x-u1(x,y)dt,y-u2(x,y)dt).
Remark: Note that convect is a nonlocal operator. The statement f=convect(f,u,v,dt) would give an incorrect result because it modifies f at some vertex where the old value is still needed later. It is necessary to do
g=convect(f,u,v,dt); f=g;
Its purpose is to define a boundary condition for a Partial Differential Equation (PDE).
The general syntax is
onbdy(ib1, ib2,...) id(u)+<expression>*dnu(u) = g onbdy(ib1, ib2,...) u = g
where ib's are boundary identification numbers, <expression> is a generic expression and g a generic function. The term id(u) may be absent as in -dnu(u)=g. dnu(u) represents the conormal of the PDE, i.e. Cgrad(u).n when the PDE operator is a*u-div(C.grad(u)).
The syntax for pde is
pde(u) Cop1(u)[*/]exp1 C op2(u)[*/]exp2...=exp3
It defines a partial differential equation with non constant coefficients where op is one of the operator:
and where [*/] means either a * or a / and similarly for C. Note that the expressions are necessarily AFTER the operator while in practice they are between the 2 differentiations for laplace...dyy. Thus laplace(u)*(x+y) means div((x+y)grad(u)).Similarly dxy(u)*f means dx(f*dy(u)).
The syntax for a single unknown function u solution of a PDE is
solve(u) begin onbdy()...; onbdy()...; ...; pde(u)... end;
For 2-systems and the use of solve(u,v), see the section 2-Systems. It defines a PDE and its boundary conditions. It will be solved by the Finite Element Method of degree 1 on triangles and a Gauss factorization.
Once the matrix is built and factorized solve may be called again by solve(u,-1)...; then the matrix is not rebuilt nor factorized and only a solution of the linear system is performed by an up and a down sweep in the Gauss algorithm only. This saves a lot of cpu time whenever possible. Several matrices can be stored and used simultaneously, in which case the sequence is
solve(u,i)...; ... solve(u,-i)...;
where i is a scalar variable (not an array function).
However matrices must be constructed in the natural order: i=1 first then i=2.... after they can be re-used in any order. One can also re-use an old matrix with a new definition, as in
solve(u,i)...; ... solve(u,i)...; solve(u,ħi)...;
Notice that solve(u) is equivalent to solve(u,1).
Remark: 2-Systems have their own matrices, so they do not count in the previous ordering.
Before going to systems make sure that your 2 pde's are indeed coupled and that no simple iteration loop will solve it, because 2-systems are significantly more computer intensive than scalar equations.
Systems with 2 unknowns can be solved by
solve(u,v) begin onbdy(..) ...dnu(u)...=.../* defines a bdy condition for u */ onbdy(..) u =... /* defines a bdy conditions for v */ pde(u) ... /* defines PDE for u */ onbdy(..)<v=... or ...dnu(v)...> /* defines bdy conditions for v */ pde(v) ... /* defines PDE for u */ end;
The syntax for solve is the same as for scalar PDEs; so solve(u,v,1) is ok for instance. The equations above can be in any orders; several onbdy() can be used in case of multiple boundaries...
The syntax for onbdy is the same as in the scalar case; either Dirichlet or mixed-Neumann, but the later can have more than one id() and only one dnu().
Dirichlet is treated as if it was mixed Neumann with a small coefficient. For instance u=2 is replaced by dnu(u)+1.e10*u=2.e10, with quadrature at the vertices. Conditions coupling u,v are allowed for mixed Neumann only, such as id(u)+id(v)+dnu(v)=1. (As said later this is an equation for v).
In solve(u,v,i) begin .. end; when i>0 the linear system is built factorized and solved. When i<0, it is only solved; this is useful when only the right hand side in the boundary conditions or in the equations have change. When i<0, i refers to a linear system i>0 of SAME TYPE, so that scalar systems and 2-systems have their own count.
Remark: saveall('filename',u,v) works also.
The syntax for pde() is the same as for the scalar case. Deciding which equation is an equation for u or v is important in the resolution of the linear system (which is done by Gauss block factorization) because some block may not be definite matrices if the equations are not well chosen.
Corners where boundary conditions change from Neumann to Dirichlet are ambiguous because Dirichlet conditions are assigned to vertices while Neumann conditions should be assigned to boundary edges; yet Gfem does not give access to edge numbers. Understanding how these are implemented helps overcome the difficulty.
All boundary conditions are converted to mixed Fourier/Robin conditions:
id(u) a + dnu(u) b = c;
For Dirichlet conditions a is set to 1.0e12 and c is multipled by the same; for Neumann a=0. Thus Neumann condition is present even when there is Dirichlet but the later overrules the former because of the large penalty number. Functions a,b,c are piecewise linear continuous, or discontinuous if precise is set.
In case of Dirichlet-Neumann corner (with Dirichlet on one side and Neumann on the other) it is usually better to put a Dirichlet logic at the corner. But if fine precision is needed then the option precise can guarrantee that the integral on the edge near the corner on the Neumann side is properly taken into account because then the corner has a Dirichlet value and a Neumann value by the fact that functions are discontinuous.
Within a program the keyword plot(u) will display u.
Instruction save('filename',u) will save the data array u on disk. If u is a scalar variable then the (single) value of u is appended to the file (this is useful for time dependent problems or any problem with iteration loop.).
Instruction savemesh('filename') will save the triangulation on disk
Similarly for reading data with load('filename',u) and loadmesh('filename').
The file must be in the default directory, else it won't be found. The file format is best seen by opening them with a text editor. For a data array f it is:
ns f[0] .... f[ns-1]
(ns is the number of vertices)
If f is a constant, its single value is appended to the end of the file; this is useful for time dependant problems or any problem with iteration loop.
If precise is set still the function stored by save is interpolated on the vertices as the P^1 continuous function given by mass lumping (see above).
For triangulations the file format is (nt = number of triangles):
ns nt q[0].x q[0].y ib[i] ... q[n-1].x q[n-1].y ib[n-1] me[0][0] me[0][1] me[0][2] ngt[0] ... me[n-1][0] me[n-1][1] me[n-1][2] ngt[n-1]
Remark: Gfem uses the Fortran standard for me[][] and numbers the vertices starting from number 1 instead of 0 as in the C-standard. Thus in C-programs one must use me[][]-1.
Remark: other formats are also recognized by freefem via their file name extentions for our own internal use we have defined .amdba and .am_fmt. You can do the same if your format is not ours.
loadmesh('mesh.amdba'); /* amdba format (Dassault aviation) */
loadmesh('mesh.am_fmt'); /* am_fmt format of MODULEF */
The purpose is to solve all the data for a PDE or a 2-system with only one instruction. It is meant for those who want to write their own solvers.
The syntax is:
saveall('file_name', var_name1,...)
The syntax is exactly the same as that of solve(,) except that the first parameter is the file name. The other parameters are used only to indicate to the interpreter which is/are the unknown function.
The file format for the scalar eq (laplace is decomposed on nuxx,nuyy)
u=p if Dirichlet c u+dnu(u)=g if Neumann b u-dx(nuxx dx(u))-dx(nuxy dy(u))-dy(nuyx dx))-dy(nuyy dy(u)) + a1 dx(u) + a2 dy(u) =f
is that each line has all the values for x,y being a vertex: f, g, p, b, c, a1, a2, nuxx, nuxy, nuyx, nuyy.
The actual routine is in C++
int saveparam(fcts *param, triangulation* t, char *filename, int N)
{
int k, ns = t->np;
ofstream file(filename);
file<<ns<<" "<<N<<endl;
for(k=0; k<ns; k++)
{
file << (param)->f[k]<<" " ; file<<" ";
file << (param)->g[k]<<" " ; file<<" ";
file << (param)->p[k]<<" " ; file<<" ";
file << (param)->b[k]<<" " ; file<<" ";
file << (param)->c[k]<<" " ; file<<" ";
file << (param)->a1[k]<<" " ; file<<" ";
file << (param)->a2[k]<<" " ; file<<" ";
file << (param)->nuxx[k]<<" " ; file<<" ";
file << (param)->nuxy[k]<<" " ; file<<" ";
file << (param)->nuyx[k]<<" " ; file<<" ";
file << (param)->nuyy[k]<<" " ; file<<" ";
file << endl;
}
}
The same function is used for complex coefficients, by overloading the operator <<:
friend ostream& operator<<(ostream& os, const complex& a)
{
os<<a.re<<" " << a.im<<" ";
return os;
}
For 2-systems also the same is used with
ostream& operator<<(ostream& os, cvect& a)
{
for(int i=0; i<N;i++)
os<<a[i]<<" ";
return os;
}
ostream& operator<<(ostream& os, cmat& a)
{
for(int i=0; i<N;i++)
for(int j=0; j<N;j++)
os<<a(i,j)<<" ";
return os;
}
where N=2.
A Dirichlet condition is applied whenever p[k](?). ( Dirichlet conditions with value 0 are changed to value 1e-10)
Gfem supports other interesting features:
The syntax is:
iter(j){....}
where j refers to the number of loops; j can be the result of an expression (as in iter(i*k)).
Imbedded loops are not allowed.
You can use iter with the adaptation features of Gfem (see section Mesh adaptation).
Gfem can handle complex coefficients with 4 dedicated keywords:
There is purposely no conjug function but barz=Re(z)-I*Im(z) will do.
By default all graphics display the real part. To display the imaginary part do plot(Im(f)). WARNING: failure to declare complex in the program implies all computation will be done in real, even if I is used.
WARNING: usual functions like cos, sin... work with real number only even if complex is set. Thus cos(I) is like cos(Re(I)) which is 1.
The linear systems for the PDE are solved by a Gauss complex LU factorization.
The instruction a:=scal(f,g); does a = int{f(x,y) conj(g)(x,y) dxdy} on Omega
where Omega is the triangulated domain.
Whenever there is a plot command, Gfem stops to let the user see the result. By using nowait no stop will be made;
Remark under X11: If you click the right button in the window, the next time the solver will give the hand to the plotter the program will stop.
This function is only available under integrated environments.
The last function defined by the keyword plot is displayed. It can be visualized in several fashion, one of which being a one dimensional plot along any segment defined by the mouse. Selection of this menu brings causes Gfem to waits for the user input which should be the line segment on which the function is to be displayed. Thus one should press the mouse at the beginning point then drag the mouse and release the button at the final point
It is safe to click in the window after to check that the function display is correct. What is seen is a
t -> f(x(t),y(t))
plot where [x(t),y(t)]_t is the segment drawn by the user and f is the last function displayed in the plot window.
The abscissa is the distance with the beginning point.
This keyword warns Gfem that precise quadrature formula must be used. Hence array-functions are discretized as piecewise linear discontinuous functions on the triangulation. Then all integrals are computed with 3 inner Gauss points slightly inside each triangle.
This option consumes more memory (3nt instead of ns per functions, i.e. 9 times more approximately, but still it is nothing compared with the memory that a matrix of the linear system of a PDE requires) because each array-function is stored by 3 values on each triangles.
It is a good idea to use it in conjunction with convect and/or discontinuous nonlinear terms.
exec('prog_name') will launch the application prog_name. It is useful to execute an external PDE solver for instance especially under Unix. It is not implemented for Macintosh because there is no simple way to return to MacGfem after prog_name has ended. The same can be achieved manually by a suitable combination of saveall, wait and load and simultaneous execution under multifinder.
user(what,f) calls the C++ function in a j-loop:
for (int j=0; j<nquad, j++) creal gfemuser( creal what, creal* f, int j)
Remark: An example of such gfemuser function is in fem.C; if you wish to put your own you must compile and link it. Under Unix, it is easy. Under the Macintosh system, either you use freefem, which is Gfem's kernel, or you must ask us a library version of MacGfem.
Suppose one wants to add an instruction to freefem, here is what must be done:
Here is the very simple structure we used for the nodes of the tree:
typedef struct noeud
{
Symbol symb;
creal value;
ident *name;
long junk;
char *path;
struct noeud *l1, *l2, *l3, *l4;
} noeud, *arbre;
In Gfem you can use mesh adaptation to improve the mesh qualities according to the solutions(one or more) of a first computing; it is easy to use this feature of Gfem, all you have to do is to use the reserved word adaptmesh.
This code was developped by Manollo J. Castro. It was ported for Gfem by Prud'homme christophe.
In order to simplify the documentation and to go straigth to direct applications, let's look at one example: the step. Here is the Gfem script:
nowait;
border(1,0,1,6,5) begin x:=0; y:=1-t end;
border(2,0,1,15,2) begin x:=2*t; y:=0 end;
border(2,0,1,10,2) begin x:=2; y:=-t end;
border(2,0,1,20,2) begin x:=2+3*t; y:=-1 end;
border(2,0,1,35,2) begin x:=5+15*t; y:=-1 end;
border(4,0,1,10,2) begin x:=20; y:=-1+2*t end;
border(5,0,1,35,5) begin x:=5+15*(1-t); y:=1 end;
border(5,0,1,40,5) begin x:=5*(1-t);y:=1 end;
buildmesh(1500);
nu := 0.001; dt := 0.4;
/* initial values */
u = y*(1-y)*one(y>0);
v = 0;
p = 2*nu*x*one(y>0);
un = u; vn = v;
i:=1; j:=2; k:=3;count:=1;
iter(500)
begin
if (count==100) then { /* adapt every 100 iterations */
adaptmesh(u,p);
i:=1; j:=2;k:=3;
wait;
plot(u);plot(v);plot(p);
nowait;
count:=0; /* reinit the counter */
};
/* increment the counter */
count:=count+1;
f=convect(un,u,v,dt); g=convect(vn,u,v,dt);
solve(u,i){ /*Horizontal velocity*/
onbdy(1) u = y*(1-y);
onbdy(2,5) u = 0;
onbdy(4)dnu(u)=0;
pde(u) id(u)/dt-laplace(u)*nu = f/dt -dx(p);
};
solve(v,j){ /* Vertical velocity */
onbdy(1,2,4,5) v = 0;
pde(v) id(v)/dt-laplace(v)*nu = g/dt -dy(p);
};
solve(p,k) { /* Pressure */
onbdy(1,2,5) dnu(p) = 0;
onbdy(4) p=0;
pde(p) -laplace(p)= -(dx(f) + dy(g))/dt;
};
un = u; vn = v; i:=-1; j:=-2; k:=-3;
end ;
wait;
plot(u);plot(v);plot(p);
As you can see in the example, we are doing 500 iterations and we adapt the mesh every 100 iterations. This is done with the counter count which is initialized to 1 before the loop, then it is incremented every iterations until 100. Then we use the if-then construction to test whether the counter is equal to 100 or not. If it is the case, we adapt the mesh according to two functions u,p; we put i,j,k to their positive value so that the matrices will be reconstructed and finally we reinitialize the counter to 1.
After the adpatation, the mesh is automatically redrawn(there is no need to use a trick like
w=0; plot(w);
to visualize the newly adapted mesh.
Note the new argument for border. This last argument is the reference of the first point of this border(it is optional); this allows to define precisely the references of each edges of the boundary. It is strongly advised to use this when you do mesh adaptation, otherwise you will have undefined edges and vertices references.
As you have certainly seen freefem is not user-friendly, so if you want to change certain options of this utility, you will haver to change the code and recompile it :-((.
However in gfem this tool is competely integrated via a dialog box. All you have to do is to click on the menu Execute, then choose the Adaptation mesh setup. A dialog box will be popped up. The options are selected by default like in freefem.
You can choose between isotropic or anisotropic meshes. One will be generate meshes with no main stretching direction, the other will have one so that the triangle will be stretched along the stream.
The mesh tool offers four kind of interpolation
They are parameters which are not often changed, their value by default are good. However you might want to change them as long as you know what you are doing.
You can type the minimal angle to determine where are the corners; There are others options, but it is better to leave these untouched.
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);
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 };
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);
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 */
/* 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);
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
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;
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);
/* 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);
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);
/* 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));