File: generalities.tex

package info (click to toggle)
freefem 3.5.8-4.2
  • links: PTS
  • area: main
  • in suites: squeeze
  • size: 2,476 kB
  • ctags: 930
  • sloc: sh: 10,708; cpp: 10,323; perl: 121; makefile: 115
file content (1355 lines) | stat: -rw-r--r-- 45,241 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
%
% $Id: generalities.tex,v 1.4 2001/10/20 23:57:17 prudhomm Exp $
%
% SUMMARY:      
% USAGE:        
%
% AUTHOR:       Christophe Prud'homme <prudhomm@mit.edu>
% 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(<max\_number\_of\_vertices>, <flag>)}\\
  where \texttt{<flag>} =
  
    $$    \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)+<expression>*dnu(u) = g
onbdy(ib1, ib2,...) u = g
\end{Verbatim}
where \verb+ib+'s are boundary identification numbers, \verb+<expression>+  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<sect>  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(..)<v=... or ...dnu(v)...> /* 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(<unknown function list>;<test function list>
          ,<<int>>) <<instruction>> : <expression>>;
\end{Verbatim}

where
\begin{itemize}
\item 
\verb+<unknown function list>+ and
\item 
\verb+<test function list>+
are one or two function names separated by a comma.
\item
\verb+<int>+ is a positive or negative integer
\item
instruction is one instruction or more if they are
   enclosed within begin end or \texttt{\{\}}
\item
\verb+<expression>+ 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(<filename>, <function name> 
    [,<variable counter: integer or converted to integer>])
load(<filename>, <function name> 
    [,<variable counter: integer or converted to integer>])
savemesh(<filename>[,<variable counter: 
                      integer or converted to integer>])
loadmesh(<filename>[,<variable counter: 
                      integer or converted to integer>])
\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{'<generic filename>-<counter>'}.
%% 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<<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;
  }
}
\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<<a.re<<" " << a.im<<"             ";  
   return os;   
 }
\end{Verbatim}
For 2-systems also the same is used with
\begin{Verbatim}
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;   
  }
\end{Verbatim}         
where \verb+N=2+ .           

A Dirichlet condition is applied whenever p[k](?).  ( Dirichlet
conditions with value \verb+0+ are changed to value \verb+1e-10+ )


\section{Other features}

\textbf{Gfem} supports other interesting features:

\subsection{Iter}

\index{iter@\texttt{iter}}
The syntax is: 
\begin{Verbatim}
iter(j){....}
\end{Verbatim}
where j refers to the number of loops; j can be the result of an expression
(as in \texttt{iter(i*k)}).

Imbedded loops are not allowed. You can use \verb+iter+  with the adaptation
features of \textbf{Gfem}.


%% node  Complex numbers, Scal(), Iter, Other features
\subsection{Complex numbers}
%%  node-name,  next,  previous,  up

\index{complex@\texttt{complex}}
\textbf{Gfem} can handle complex coefficients with  4 dedicated keywords:
\begin{itemize}
\item
\verb+complex+ : to tell Gfem that complex number will be used.
When it is used it must be located at the beginning of the program 
before any  function declarations, otherwise the results will be incorrect.
It can appear more than once in the program but only the first occurrence
counts.

\item \texttt{I}: for \texttt{sqrt(-1)} .

\item \texttt{Re} : for the real part.

\item \texttt{Im} : for the imaginary part.
\end{itemize} 

There is purposely no conjug function but \texttt{barz=Re(z)-I*Im(z)}  will
do.

By default all graphics display the real part. To display the imaginary
part do \texttt{plot(Im(f))}.
 
The functions implemented for complex numbers are:
\begin{itemize}
\item $cos$
\item $sin$
\item $z^x$ where $z$ is complex and $x$ is a float
\end{itemize}

The linear systems for the PDE are solved by a Gauss complex LU 
factorization.

\textbf{WARNING}: failure to declare \verb+complex+  in the program implies all
computation will be done in real, even if \verb+I+  is used. 

%% node  Scal(), Wait, Complex numbers, Other features
\subsection{Scal()}
%%  node-name,  next,  previous,  up

\index{scal@\texttt{scal}}
The instruction \verb+a:=scal(f,g);+  does

  $$ a = \int_\Omega f(x,y) \overline{g}(x,y) dxdy  $$
where $\Omega$ is the triangulated domain.

%% node  Wait, One Dimensional Plots, Scal(), Other features
\subsection{Wait, nowait, changewait}
%%  node-name,  next,  previous,  up

Whenever there is a \verb+plot+ command, \textbf{Gfem}  stops to let the user see
the result. By using \verb+nowait+  no stop will be made; 
\begin{itemize}
\item 
\texttt{wait}  turns back the stop option on.
\item
\texttt{changewait} toggles the option from on to off or off to on.
\end{itemize}

\textbf{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.

%% node  One Dimensional Plots, Precise, Wait, Other features
\subsection{One Dimensional Plots}
%%  node-name,  next,  previous,  up

This function is \textbf{only available} under integrated environments.

The last function defined by the keyword \verb+plot+  is displayed.  It can
be visualized in several fashion, one of which being  a \textbf{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  \textbf{
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 \to f(x(t),y(t))$$

Plot where $[x(t),y(t)]_t$ is the segment drawn by the user and \verb+f+  is
the last function displayed in the plot window.  

The abscissa is the distance with the beginning point.

%% node   Precise, Exec(), One Dimensional Plots, Other features
\subsection{Precise}
%%  node-name,  next,  previous,  up

\index{precise@\texttt{precise}}
This keyword warns \textbf{Gfem} that precise quadrature formula must be
used. Hence array-functions are discretized as piecewise linear
\textbf{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 \verb+convect+
\index{convect@\texttt{convect}} and/or discontinuous nonlinear terms.

%% node  Exec(),  , Precise, Other features
\subsection{Exec(), user(), how to link an external function to Gfem}
%%  node-name,  next,  previous,  up

\index{exec@\texttt{exec}}
\index{user@\texttt{user}}

\verb+exec('prog_name')+ will launch the application \verb+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 \texttt{progname} has ended.  
The same can be achieved  manually by a suitable combination of \verb+saveall+, \verb+wait+  and
\verb+load+  and simultaneous execution under multifinder.

\verb+user(what,f)+  calls the C++ function in a j-loop:
\begin{Verbatim}
for (int j=0; j<nquad, j++)
creal gfemuser( creal what, creal* f, int j)
\end{Verbatim}
\begin{itemize}

\item
\verb+creal+ is a scalar (\verb+float+ ) or a complex number if \verb+complex+  has
been set; \verb+nquad=3*nt+ if \verb+precise+  is set and \verb+ns+  otherwise. 

\item
\verb+what+  is intended for users who need several such functions.  
Then all can be put in the super function \verb+user+  and selection is by
an if statement on \verb+what+ .

\item
Within \verb+gfemuser+  access to all global variables are of course possible:
the triangulation (\verb+ns,nt, me, q, ng,ngt, area...+ ) ... refer to the file 
\verb+fem.C+  for more details.

\end{itemize}

\textbf{Remark}: An example of such gfemuser function is in \verb+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.

%% language internals
%% node  Language internals,  , Other features, Generalities
\section{Language internals}
%%  node-name,  next,  previous,  up

\begin{itemize}
\item
Like most interpreters it has a lexical and a syntaxic part.  
In the lexical part the source is broken into tokens and recognized as 
symbols (see the enum symbol in the source file \texttt{lexical.h}).
For instance if the first character is a digit then it is a number and the
symbol type associated is \texttt{cste}.  This job is done by  function
\texttt{nextsym}.In addition it constructs a table of constants and variables
(which for convenience contains also the reserved words of the language).
\item
The lexical analyzer is a function called by the syntax analyzer. 
Hence the second is the main routine in the program except for a few 
initialization; it's name is \verb+instruction+ .The syntax analysis is
driven by the syntax rules because the language is LL(1).  Thus there
is C-function for each non-terminal.
\item
The program does not generate an object code but a tree.
For example,  parsing  \verb+x * y+ would generate a tree with root \verb+ * +   
and two branches \verb+x+ and \verb+y+ . Trees are C-struct with four pointers
for the 4 branches (here two would point to NULL) and a symbol for the
operator. The C-function which builds  trees is called \verb+plante+ .          
In the end the program is transformed into a full tree and 
to execute the program there is only one thing to do: evaluate
the operator at the root of the tree.
\item
The evaluation of the program is done by the C-function \verb+eval+ . It
looks at the root symbol and perform the corresponding operation. Here
it would do: \emph{return "value pointed by L1"*"value pointed by L2" if L1
and L2 where the addresses of the two branches}. 
\item
The art of the game is to associate a tree to each operation. For 
example when  the value of the variable \verb+x+  is required, this is also done 
by a tree which has   the operator "oldvar" as root.  The trickiest of all is 
the compound instruction \verb+{...;....};+ .  
Here \verb+{+  is considered as an operator with one branch on the current 
instruction and one branch on the next one.  Similarly for the 
\verb+ if...then... else+  instruction.

\end{itemize}

Suppose one wants to add an instruction to \textbf{freefem}, here is what must be
done:
\begin{itemize}
\item
Make sure the syntax is LL(1) and does not conflict with the old one. 
\item
Add reserved words to the table with \verb+installe+ . As there will be a
new Symbol, update the list of symbols (an enum structure). Add the
C-functions for the syntax analysis according  to the diagrams. Modify
eval by adding to the \verb+switch+  the new case. 
\end{itemize}

Here is the very simple structure we used for the nodes of the tree:
\begin{Verbatim}
typedef struct noeud
{
  Symbol symb; 
  creal value;
  ident *name;
  long  junk;
  char  *path;
  struct noeud *l1, *l2, *l3, *l4;
} noeud, *arbre;
\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:
%