## File: generalities.tex

package info (click to toggle)
freefem 3.5.8-4.2
 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355 % % $Id: generalities.tex,v 1.4 2001/10/20 23:57:17 prudhomm Exp$ % % SUMMARY: % USAGE: % % AUTHOR: Christophe Prud'homme % ORG: MIT % E-MAIL: prudhomm@mit.edu % % ORIG-DATE: 8-Feb-97 at 16:45:11 % LAST-MOD: 20-Oct-01 at 15:45:01 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{The Gory Details} \label{cha:general} \section{Programs} \label{sec:programs} \textbf{Gfem} is 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 \texttt{begin} can be replaced by \texttt{\{} and \texttt{end} by \texttt{\}}. \textbf{C programmers}: caution the syntax is "...\};" while most C constructs use " ...;\}" \textbf{ Example 1}: Triangulates a circle and plot $f = x*y$ \begin{Verbatim} border(1,0,6.28,20) begin x:=cos(t); y:=sin(t); end; buildmesh(200); f=x*y; plot(f); \end{Verbatim} \textbf{ Example 2}: on the circle solve the Dirichlet problem -$\Delta(u) = x*y$ with $u=0$ on $\partial \Omega$ \begin{Verbatim} 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); \end{Verbatim} \section{List of reserved words} \begin{table}[htbp] \begin{center} \begin{tabular}[c]{|l|l|} \hline \multicolumn{1}{|c|}{Keywords} & \multicolumn{1}{|c|}{Explanations} \\ \hline begin, \{ & Begin a new block \\ \hline end, \} & End a block \\ \hline if, then, else, or, and,iter & Conditionnals and Loops \\ \hline x,y,t,ib,iv,region,nx,ny & Mesh variables \\ \hline log, exp, sqrt, abs, min, max & \\ sin, asin, tan, atan, cos, acos & Mathematical Functions\\ cosh, sinh, tanh & \\ \hline I, Re, Im & Complex Numbers \\ complex & Enter Complex Number Mode \\ \hline buildmesh, savemesh, loadmesh, adaptmesh & Mesh related Functions \\ & build, save , load and mesh adaptation\\ \hline one, scal, dx, dy, dxx, dyy, dxy, dyx, convect & Mathematical operators \\ & can be called wherever you want \\ \hline solve & Enter Solver Mode \\ pde, id, dnu laplace, div onbdy & Solver Functions\\ plot, plot3d & graphical functions: plot isolines in 2D\\ & and Elevated Surface in 3D \\ \hline save, load, saveall & Saving and Loading Functions \\ \hline & Change the Wait State: if \texttt{wait}\\ wait, nowait, changewait & then the user must click in \\ & the window to continue \\ \hline precise, halt, include, evalfct, exec, user & Miscellaneous Functions \\ \hline \end{tabular} \caption{Gfem Keywords} \label{tab:1} \end{center} \end{table} \section{Building a mesh} \subsection{Triangulations} To create a triangulation you must either \begin{itemize} \item Open a project \item Read an old triangulation stored in text format \item Execute a program which contains the keyword \verb+buildmesh+ \item Create one by hand drawing the boundary of $\Omega$ and activate the menu \verb+Triangulate+ (Macintosh only). \end{itemize} In integrated environments, once created, triangulations can be displayed, stored on disk in \textbf{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). \textbf{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. %% node Border(), Geometric variables, Triangulations, Building a mesh \subsection{Border(), buildmesh(), polygon()} %% node-name, next, previous, up \index{border@\texttt{border}} \index{buildmesh@\texttt{buildmesh}} \index{polygon@\texttt{polygon}} Use it to triangulate domain defined by its boundary. The syntax is \begin{Verbatim} border(ib,t_min,t_max,nb_t) begin ...x:=f(t); ...y:=g(t)... end; buildmesh(nb_max); \end{Verbatim} where each line with \verb+border+ could be replaced by a line with \verb+polygon+ \begin{Verbatim} polygon(ib,'file_name'[,nb_t]); \end{Verbatim} where \verb+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 \verb+t+ and the two coordinates must be \verb+x+ and \verb+y+ . When the parameter goes from \verb+t_min+ to \verb+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 \verb+border+ , and have inner dividing curves; \verb+nb_t+ points are put on the boundary with values $t = t_{min} + i * (t_{max} - t_{min}) / (nb_t-1)$ where \verb+i+ takes integer values from \verb+0+ to \verb+nb_t-1+ . The triangulation is created by a Delaunay-Voronoi method with \verb+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 \verb+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 \verb+ib=0+ . The usual \verb+if... then ... else+ statement can be used with the compound statement: \verb+begin...end+ . This allows piecewise parametric definitions of complicated or polygonal boundaries. The boundary identification number \verb+ib+ can be overwritten.For example: \begin{Verbatim} 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); \end{Verbatim} Recall that \verb+begin+ and \verb+{+ is the same and so is \verb+end+ and \verb+}+. Here one side of the unit square has ib=1. The 3 other sides have ib=2. The keyword \verb+polygon+ causes a sequence of boundary points to be read from the file \verb+file_name+ which must be in the same directory 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 \begin{Verbatim} x[0] y[0] x[1] y[1] .... \end{Verbatim} each being separated by a tab or a carriage return. The last parameter \verb+nb_t+ is optional; it means that each segment will be divided into \verb+nb_t+1+ equal segments (i.e. \verb+nb_t+ points are added on each segments). For example \begin{Verbatim} polygon(1,'mypoints.dta',2); buildmesh(100); \end{Verbatim} with the file mypoints.dta containing \begin{Verbatim} 0. 0. 1. 0. 1. 1. 0. 1. 0. 0. \end{Verbatim} triangulates the unit square with 4 points on each side and gives \texttt{ib=1} to its boundary. Note that \texttt{polygon(1,'mypoints.dta')} is like \texttt{polygon(1,'mypoints.dta',0)}\index{polygon@\texttt{polygon}}. \subsubsection{\texttt{buildmesh} and domain decomposition} There is a problem with \texttt{buildmesh}\index{buildmesh@\texttt{buildmesh}} when doing domain decomposition: by default \textbf{Gfem} swap the diagonals\index{diagonal swaping} at the corners of the domain if the triangle has two boundary edges. This will lead to bad domain decomposition\index{domain decomposition} at the sub-domain interfaces. \par To solve this, there is a new flag for \texttt{buildmesh} which is optional: \begin{center} \texttt{buildmesh(, )}\\ where \texttt{} = $$\left\{ \begin{array}[l]{l} = 0 \mbox{ classic way: do diagonal swaping}\\ = 1 \mbox{ domain decomposition: no diagonal swaping} \end{array} \right.$$ \end{center} %% node Geometric variables, Regions, Border(), Building a mesh \subsection{Geometric variables, inner and region bdy, normal} %% node-name, next, previous, up \begin{itemize} \item \verb+x,y+ refers to the coordinates in the domain \item \verb+ib+ refers to the boundary identification number; it is zero inside the domain. \item \verb+nx+ and \verb+ny+ refer to the x-y components of the normal on the boundary vertices; it is zero on all inner vertices. \item \verb+region+ refers to the domain identification number which is itself based on an internal number, ngt, assigned to each triangle by the triangulation constructor. \end{itemize} 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. %% node Regions, , Geometric variables, Building a mesh \subsection{Regions} %% node-name, next, previous, up A sub-domain is a piece of the domain which is simply connected and delimited by boundaries. Each sub-domain has a \textbf{region} number assigned to it. This is done by \textbf{Gfem}, not by the user. Every time \verb+border+ is called, an internal number \verb+ngt+ is incremented by 1. Then when the key word \verb+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 \verb+ngt+. The number \textbf{region} is a piecewise linear continuous interpolation at the vertices of \verb+ngt+. To be exact, the value of \textbf{region} at a vertex $(x_0,y_0)$ is the value of \texttt{ngt} at $(x_0,y_0-10^{-6})$, except if \texttt{precise} is set in which case \textbf{region} is equal to \texttt{ngt}. \section{Functions} \subsection{Functions and scalars} Functions are either read or created. \begin{itemize} \item Functions can be read from a file if its values at the vertices of the triangulation are stored in text format. (Open a .dta example with a text editor to see the format). \item Functions can be created by executing a program. An instruction like \verb+f=x*y+ really means that $f(x,y)=x*y$ for all \verb+x+ and \verb+y+. Here \verb+x+ and \verb+y+ refer to the coordinates in the domain represented by the triangulation. \item Functions can be created with other previously defined functions such as in \verb+g=sin(x*y); f=exp(g);+ . \item Four other variables can be used besides \verb+x, y, iv, t+ : \verb+ nx, ny, ib, region+. \end{itemize} Most usual functions can be used: \begin{Verbatim} max, min, abs, atan, sqrt, cos, sin, tan, acos, asin, one, cosh, sinh, tanh, log, exp \end{Verbatim} \verb+one(x+y<0)+ for instance means \verb+1+ if \verb+x+y<0+ and \texttt{0} otherwise. Operators: \begin{Verbatim} and, or, < , <=, < , >=, ==, +, -, *, /, ^ x^2 means x*x \end{Verbatim} Functions created by a program are displayed only if the key word \verb+plot()+ or \verb+plot3d()+ is used ( here \verb+plot(f)+ ). Derivatives of functions can be created by the keywords \verb+dx()+ and \verb+dy()+ . Unless \verb+precise+ is set, they are interpolated so the results is also continuous piecewise linear (or discontinuous when precise is set). Similarly the convection operator \texttt{convect(f,u1,u2,dt)} \index{convect@\texttt{convect}} 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 \verb+:=+ is useful to create them, as in \begin{Verbatim} a:= (1+sqrt(5))/2; f= x*cos(a*pi*y); \end{Verbatim} Here \verb+f+ is a function, \verb+a+ is a scalar and \verb+pi+ is a (predefined) a scalar. It is possible to evaluate a function at a point as in \verb+a:=f(1,0)+ Here the value of \verb+a+ will be \verb+1+ because \verb+f(1,0)+ means \verb+f+ at \verb+x=1+ and \verb+y=0+ . %% node Building functions, Value of a function at one point, Functions and scalars, Functions \subsection{Building functions} %% node-name, next, previous, up There are 6 predefined functions: \verb+x,y,ib,region, nx, ny+ . \begin{itemize} \item The coordinate \verb+x+ is horizontal and \verb+y+ is vertical. \item $\mbox{\texttt{ib}} = \left\{ \begin{array}[l]{l} 0 \mbox{ inside } \Omega \\ > 0 \mbox{ on } \partial \Omega \end{array} \right.$ On $\partial \Omega$ it is equal to the boundary identification number. \end{itemize} The usual \texttt{if... then ... else} statement can be used with an important restriction on the logical expression which must return a \textbf{scalar} value: \begin{Verbatim} if( logical expression) then { statement; ....; statement; } else { ..... }; \end{Verbatim} The \texttt{logical expression} controls the \texttt{if} by its return being 0 or >0, it is evaluated only once (i.e. with \texttt{x}, \texttt{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 \verb+:=+ tells the compiler not to allocate a data array to this variable. Thus \verb+v=sin(a*pi*x);+ generates an array for \verb+v+ but no array is assigned to \verb+a+ in the statement \verb+a:=2+ . \subsection{Value of a function at one point} If \verb+f+ has been defined earlier then it is possible to write \verb+a:=f(1.2,3.1);+ Then a has the value of \verb+f+ at \verb+x=1.2+ and \verb+y=3.1+ . It is also allowed to do \begin{Verbatim} x1:=1.2; y1:=1.6; a:=f(x1,2*y1); save('a.dta',a); \end{Verbatim} \textbf{Remark:} Recall that ,\verb+a+ being a scalar, its value is \textbf{appended} to the file \texttt{a.dta}. %% node Special functions \subsection{Special functions:dx(), dy(), convect()} \index{convect@\texttt{convect}} \index{dx@\texttt{dx}} \index{dy@\texttt{dy}} \verb+dx(f)+ is the partial derivative of \verb+f+ with respect to \verb+x+ ; the result is piecewise constant when \verb+precise+ is set and interpolated with mass lumping as a the piecewise linear function when precise is not set. %$$dx(u)^i = {\frac{1}{\Sigma_i\Sigma_{k,q^i in T_k} u^k |T_k|}}$$ %where $-q^i$ is a vertex,$dx(u)^i$ the value at this vertex, $\Sigma_k$ %is the area of all the triangles which have $q^i$ for vertex. $-T_k$ is %a triangle, $|T_k|$ its area, $u^k$ the mean of u on this triangle. Note that \verb+dx()+ is a non local operator so statements like \verb+f=dx(f)+ would give the wrong answer because the new value for \verb+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; \verb+convect+ \index{convect@\texttt{convect}} provides a tool for Lagrangian upwinding. By \verb+g=convect(f,u,v,dt)+ \textbf{Gfem} construct an approximation of $$f(x-u1(x,y)dt,y-u2(x,y)dt)$$ Recall that when $$\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} + v\frac{\partial f}{\partial y} = lim_{dt \to 0}\frac{f(x,y,t) - f(x-u(x,y)dt,y-v(x,y)dt,t-dt)}{dt}$$ Thus to solve $$\frac{\partial f}{\partial t} + u\frac{\partial f}{\partial x} + v\frac{\partial f}{\partial y} -div(\mu grad f) = g,$$ in a much more stable way that if the operator \verb+dx(.)+ and \verb+dy(.)+ were use, the following scheme can be used: \begin{Verbatim} iter(n) begin id(f)/dt - laplace(f)*mu =g + convect(oldf,u,v,dt)/dt; oldf = f end; \end{Verbatim} \textbf{Remark:} Note that \verb+convect+ \index{convect@\texttt{convect}} is a nonlocal operator. The statement \texttt{f = convect(f,u,v,dt)} would give an incorrect result because it modifies \verb+f+ at some vertex where the old value is still needed later. It is necessary to do \begin{Verbatim} g=convect(f,u,v,dt); f=g; \end{Verbatim} %% node Integrals, Solving an equation, Functions, Generalities \section{Global operators} \label{sec:globop} It is important to understand the difference between a global operator such as \verb+dx()+ and a local operator such as \verb+sin()+ . To compute \verb+dx(f)+ at vertex \verb+q+ we need \texttt{f} at all neighbors of \verb+q+. Therefore evaluation of \texttt{dx(2*f)} require the computation of \verb+2*f+ at all neighbor vertices of \verb+q+ before applying \texttt{dx()} ; but in which memory would the result be stored? Therefore Gfem does not allow this and forces the user to declare a function \verb+g =2*f+ before evaluation of \verb+dx(g)+ ; Hence in \begin{Verbatim} g = 2*f; h = dx(g) * dy(f); \end{Verbatim} the equal sign forces the evaluation of \verb+g+ at all vertices, then when the second equal signs forces the evaluation of the expression on the right of h at all vertices , everything is ready for it. Global operators are \begin{Verbatim} dx(), dy(), convect(), intt[], int()[] \end{Verbatim} Example of forbidden expressions: \begin{Verbatim} intt[f+g], dx(dy(g)), dx(sin(f)), convect(2*u...) \end{Verbatim} \section{Integrals} %% node-name, next, previous, up \begin{itemize} \item \textbf{effect}: \verb+intt+ returns a complex or real number, an integral with respect to \verb+x,y+ \verb+int+ returns a complex or real number, an integral on a curve \item \textbf{syntax}: \verb+intt[f]+ or \verb+intt(n1)[f]+ or \verb+intt(n1,n2)[f]+ or \verb+intt(n1,n2,n3)[f]+ \verb+int(n1)[f]+ or \verb+int(n1,n2)[f]+ or \verb+int(n1,n2,n3)[f]+ where \verb+n1,n2,n3+ are boundary or subdomain identification numbers and where \verb+f+ is an array function. \begin{Verbatim} border(1)... end; /* a border has number 1 */ ... buildmesh(...); f = 2 * x; /* * nx,ny are the components of the boundary normal */ g = f * (nx + ny); /* * can't do r:= int[2*x] */ r:= int[f]; s:=int(1)[g]; /* * this is the only way to display the result */ save('r.dta',r); save('s.dta',s); \end{Verbatim} \item \textbf{Restrictions}: \verb+int+ and \verb+intt+ are global operators, so the values of the integrands are needed at all vertices at once, therefore you can't put an expression for the integrand, it must be a function. Be careful to check that the region number are correct when you use \verb+intt(n)[f]+ . Unfortunately \textbf{freefem} does not store the edges numbers. Hence there are ambiguities at vertices which are at the intersections of 2 boundaries. The following convention is used: \verb+int(n)[g]+ computes the integral of \verb+g+ on all segments of the boundary (both ends have id boundary number !=0) with one vertex boundary id number = n. (Remember that you can control the boundary id number of the boundary ends by the order in which you place the corresponding \verb+border+ call or by an extra argument in \verb+border+ ) \end{itemize} %% node Solving an equation, Results, Integrals, Generalities \section{Solving an equation} %% node-name, next, previous, up %% Onbdy():: %% Pde():: %% Solve():: %% 2-Systems:: %% Boundary conditions at corners:: %% node Onbdy(), Pde(), Solving an equation, Solving an equation \subsection{Onbdy()} Its purpose is to define a boundary condition for a Partial Differential Equation (PDE). The general syntax is \begin{Verbatim} onbdy(ib1, ib2,...) id(u)+*dnu(u) = g onbdy(ib1, ib2,...) u = g \end{Verbatim} where \verb+ib+'s are boundary identification numbers, \verb++ is a generic expression and \texttt{g} a generic function. The term \verb+id(u)+ may be absent as in \verb+-dnu(u)=g+ . \verb+dnu(u)+ represents the conormal of the PDE, i.e. $$\pm\vec\nabla u . n$$ when the PDE operator is $$a * u - \vec\nabla.(\pm\vec\nabla u)$$ %% node Pde(), Solve(), Onbdy(), Solving an equation \subsection{Pde()} %% node-name, next, previous, up The syntax for pde is \begin{Verbatim} pde(u) [+-] op1(u)[*/]exp1 [+-] op2(u)[*/]exp2...=exp3 \end{Verbatim} It defines a partial differential equation with non constant coefficients where \verb+op+ is one of the operator: \begin{itemize} \item \verb+id()+ \item \verb+dx()+ \item \verb+dy()+ \item \verb+laplace()+ \item \verb+dxx()+ \item \verb+dxy()+ \item \verb+dyx()+ \item \verb+dyy()+ \end{itemize} and where \verb+[*/]+ means either a \verb+*+ or a \verb+/+ and similarly for $\pm$. Note that the expressions are necessarily \textbf{AFTER} the operator while in practice they are between the 2 differentiations for \verb+laplace...dyy+ . Thus \texttt{laplace(u)*(x+y)} means $\nabla.((x+y)\nabla u)$ .Similarly \texttt{dxy(u)*f} means $\frac{\partial f \frac{\partial u}{\partial y}}{\partial x}$. %% node Solve(), 2-Systems, Pde(), Solving an equation \subsection{Solve()} %% node-name, next, previous, up Solve() \index{solve@\texttt{solve}} The syntax for a single unknown function u solution of a PDE is \begin{Verbatim} solve(u) begin onbdy()...; onbdy()...; ...; pde(u)... end; \end{Verbatim} For 2-systems and the use of \verb+solve(u,v)+, see the section \verb+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 \verb+solve+ may be called again by \verb+solve(u,-1)...;+\index{solve@\texttt{solve}} 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 \begin{Verbatim} solve(u,i)...; ... solve(u,-i)...; \end{Verbatim} where \verb+i+ is a scalar variable (not an array function). However matrices must be constructed in the natural order: \verb+i=1+ first then \verb+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 \begin{Verbatim} solve(u,i)...; ... solve(u,i)...; solve(u,\pm i)...; \end{Verbatim} Notice that \verb+solve(u)+ is equivalent to \verb+solve(u,1)+\index{solve@\texttt{solve}} . \textbf{Remark:} 2-Systems have their own matrices, so they do not count in the previous ordering. %% node 2-Systems, Boundary conditions at corners, Solve(), Solving an equation \subsection{2-Systems} %% node-name, next, previous, up 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 \begin{Verbatim} 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(..) /* defines bdy conditions for v */ pde(v) ... /* defines PDE for u */ end; \end{Verbatim} The syntax for \verb+solve+ is the same as for scalar PDEs; so \verb+solve(u,v,1)+\index{solve@\texttt{solve}} is ok for instance. The equations above can be in any orders; several \verb+onbdy()+ can be used in case of multiple boundaries... The syntax for \verb+onbdy+ is the same as in the scalar case; either Dirichlet or mixed-Neumann, but the later can have more than one \verb+id()+ and only one \verb+dnu()+ . Dirichlet is treated as if it was mixed Neumann with a small coefficient. For instance u=2 is replaced by \texttt{dnu(u)+1.e10*u=2.e10} , with quadrature at the vertices. Conditions coupling \verb+u,v+ are allowed for mixed Neumann only, such as \texttt{id(u)+id(v)+dnu(v)=1}. (As said later this is an equation for \verb+v+ ). In \verb+solve(u,v,i) begin .. end;+ when \verb+i>0+ the linear system is built factorized and solved. When \verb+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 \verb+i<0+, \verb+i+ refers to a linear system \verb+i>0+ of \textbf{SAME TYPE}, so that scalar systems and 2-systems have their own count. \textbf{Remark:} \verb+saveall('filename',u,v)+ works also. The syntax for \verb+pde()+ is the same as for the scalar case. Deciding which equation is an equation for \verb+u+ or \verb+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. \begin{itemize} \item A boundary condition like \verb+ onbdy(...) ... dnu(u) ... = ...; + is a boundary condition associated to \verb+u+, even if it contains \verb+id(v)+ . \item Obviously a boundary condition like \verb+ onbdy(...) u...=...;+ is also associated with \verb+u+ (the syntax forbids any \verb+v+ -operator in this case). \item If \verb+u+ is the array function in a \verb+pde(u)+ then what follows is the PDE associated to \verb+u+ . \end{itemize} %% node Boundary conditions at corners, , 2-Systems, Solving an equation \subsection{Boundary conditions at corners} %% node-name, next, previous, up 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 \textbf{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: \begin{Verbatim} id(u) a + dnu(u) b = c; \end{Verbatim} For Dirichlet conditions a is set to \verb+1.0e12+ and \verb+c+ is multiplied by the same; for Neumann \verb+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 \verb+a,b,c+ are piecewise linear continuous, or discontinuous if \verb+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 \verb+precise+ can guarantee 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. %% node Results, Other features, Solving an equation, Generalities \subsection{Weak formulation} \label{sec:varform} The new keyword \texttt{varsolve} allows the user to enter PDEs in weak form. Syntax: \begin{Verbatim} varsolve(; ,<>) <> : >; \end{Verbatim} where \begin{itemize} \item \verb++ and \item \verb++ are one or two function names separated by a comma. \item \verb++ is a positive or negative integer \item instruction is one instruction or more if they are enclosed within begin end or \texttt{\{\}} \item \verb++ is an expression returning a real or complex number \end{itemize} We have used the notation \verb+<< >>+ whenever the entities can be omitted. Examples \begin{Verbatim} varsolve(u;w) /* Dirichlet problem -laplace(u) =x*y */ begin onbdy(1) u = 0; f = dx(u)*dx(w) + dy(u)*dy(w) g = x*y; end : intt[f] - intt[g,w]; varsolve(u;w,-1) /* same with prefactorized matrix */ begin onbdy(1) u = 0; f = dx(u)*dx(w) + dy(u)*dy(w) g = x*y; end : intt[f] - intt[g,w]; varsolve(u;w) /* Neuman problem u-laplace(u) = x*y */ begin f = dx(u)*dx(w) + dy(u)*dy(w) -x*y; g = x; end : intt[f] + int[u,w] - int(1)[g,w]; varsolve(u,v;w,s) /* Lame's equations */ begin onbdy(1) u=0; onbdy(1) v=0; e11 = dx(u); e22 = dy(v); e12 = 0.5*(dx(v)+dy(u)); e21 = e12; dive = e11 + e22; s11w=2*(lambda*dive+2*mu*e11)*dx(w); s22s=2*(lambda*dive+2*mu+e22)*dy(s); s12s = 2*mu*e12*(dy(w)+dx(s)); s21w = s12s; a = s11w+s22s+s12s+s21w +0.1*s; end : intt[a]; \end{Verbatim} \textbf{How does it works} The interpreter evaluates the expression after the ":" for each triangle and for each 3 vertices; if there is an instruction prior the ":" it is also evaluated similarly. Each evaluation is done with one of the unknown and one of the test functions being 1 at one vertices and zero at the 2 others. This will give an element of the contribution of the triangle to the linear system of the problem. The right hand side is constructed by having all unknowns equal to zero and one test function equal to one at one vertex. whenever integrals appear they are computed on the current triangle only. \textbf{Note} that \verb+varsolve+ takes longer than \texttt{solve} because derivatives like \verb+dx(u)+ are evaluated 9 times instead of once. \section{Results} %% node-name, next, previous, up %% plot:: %% Saveall():: %% node plot, Saveall(), Results, Results \subsection{plot, savemesh, save, load, loadmesh} %% node-name, next, previous, up Within a program the keyword \verb+plot(u)+ will display \verb+u+ . Instruction \texttt{save('filename',u)} will save the data array \verb+u+ on disk. If \verb+u+ is a scalar variable then the (single) value of \verb+u+ is appended to the file (this is useful for time dependent problems or any problem with iteration loop.). Instruction \texttt{savemesh('filename')}\index{savemesh@\texttt{savemesh}} will save the triangulation on disk. Similarly for reading data with \texttt{load('filename',u)} and \texttt{loadmesh('filename')}\index{loadmesh@\texttt{loadmesh}}. 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 \verb+f+ it is: \begin{Verbatim} ns f[0] .... f[ns-1] \end{Verbatim} (\verb+ns+ is the number of vertices) If \verb+f+ is a constant, its single value is \emph{appended} to the end of the file; this is useful for time dependent problems or any problem with iteration loop. If \verb+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 (\verb+nt+ = number of triangles): \index{Mesh formats!\texttt{gfem}} \begin{Verbatim} 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] \end{Verbatim} \paragraph{\textbf{Remark:}} \textbf{Gfem} uses the Fortran standard for \verb+me[][]+ and numbers the vertices starting from number 1 instead of 0 as in the C-standard. Thus in C-programs one must use \verb+me[][]-1+ . \paragraph{\textbf{Remark:}} Other formats are also recognized by freefem via their file name extensions for our own internal use we have defined \verb+.amdba+ and \verb+.am_fmt+. You can do the same if your format is not ours. \index{Mesh formats!\texttt{ambda}} \index{Mesh formats!\verb+am_fmt+} \begin{Verbatim} loadmesh('mesh.amdba'); /* amdba format (Dassault aviation) */ loadmesh('mesh.am_fmt'); /* am_fmt format of MODULEF */ \end{Verbatim} \paragraph{\textbf{Remark:}} There is an optional arguments for the functions \texttt{load}, \texttt{save}, \texttt{loadmesh}\index{loadmesh@\texttt{loadmesh}}, \texttt{savemesh}\index{savemesh@\texttt{savemesh}}. This is the 2nd or 3rd argument of these functions. Here are the prototypes: \begin{Verbatim} save(, [,]) load(, [,]) savemesh([,]) loadmesh([,]) \end{Verbatim} As an example see \texttt{nsstepad.pde} which use this feature to save the mesh and the solution at each adaptation of the mesh. This special feature allows you to save or load a generic filename with a counter, the final filename is built like this \texttt{'-'}. %% node Saveall(), , plot, Results \subsection{Saveall()} %% node-name, next, previous, up 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: \begin{Verbatim} saveall('file_name', var_name1,...) \end{Verbatim} The syntax is exactly the same as that of \verb+solve(,)+ \index{solve@\texttt{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 equation (\texttt{laplace} is decomposed on \texttt{nuxx, nuyy}) \begin{Verbatim} 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 \end{Verbatim} is that each line has all the values for \texttt{x,y} being a vertex: \texttt{f, g, p, b, c, a1, a2, nuxx, nuxy, nuyx, nuyy}. The actual routine is in C++ \begin{Verbatim} int saveparam(fcts *param, triangulation* t, char *filename, int N) { int k, ns = t->np; ofstream file(filename); file<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; } } \end{Verbatim} The same function is used for complex coefficients, by overloading the operator <<: \begin{Verbatim} friend ostream& operator<<(ostream& os, const complex& a) { os<