msgQ: comment out debug if/throw
[charm.git] / doc / fem / manual.tex
1 \documentclass[10pt]{article}
2 \usepackage{../pplmanual}
3 \input{../pplmanual}
4
5 \makeindex
6
7 \title{\charmpp\\ Finite Element Method (FEM) Framework\\ Manual}
8 \version{1.2}
9 \credits{
10 Initial version of \charmpp{} FEM Framework was developed
11 by Milind Bhandarkar with inputs from Timothy Hinrichs and Karthikeyan
12 Mahesh. The current version is almost completely rewritten by
13 Orion Lawlor. The most recent version is being called by the name ParFUM. 
14 ParFUM is short for Parallel Framework for Unstructured Meshes. 
15 This version has been mostly written by Nilesh Choudhury, Terry Wilmarth,
16 Sayantan Chakravorty and Issac Dooley.}
17
18 \begin{document}
19
20 \maketitle
21
22 \section{Introduction}
23
24 The Finite Element Method (FEM) approach is used in many engineering
25 applications with irregular domains, from elastic deformation problems to
26 crack propagation to fluid flow.  \charmpp{} is a free message-passing parallel
27 runtime system for machines from clusters of workstations to tightly-coupled
28 SMPs.  The \charmpp{} FEM framework allows you to write a parallel FEM program,
29 in C or Fortran 90, that closely resembles a serial version but includes
30 a few framework calls.
31
32 Using the FEM framework also allows you to take advantage of all the
33 features of \charmpp, including run-time load balancing,  performance
34 monitoring and visualization, and checkpoint/restart, with no additional
35 effort. The FEM framework also combines naturally with other \charmpp
36 frameworks built on TCHARM.
37
38 The FEM framework has been undergoing a wave of recent improvements. A choice to rename the new version ParFUM has been adopted.ParFUM is short for Parallel Framework for Unstructured Meshes. Section \ref{sec:ParFUM} describes some of the new features included in ParFUM that were not present in FEM.
39
40
41 \subsection{Philosophy}
42
43 The \charmpp{} FEM framework is designed to be flexible, in that it
44 provided a few very general operations, such as loading and partitioning 
45 a ``mesh.''  
46 In describing these operations, we draw on examples from structural analysis,
47 but in fact the same calls can be used for other applications, including
48 fluid dynamics or partial differential equations solvers, or
49 even general-purpose graph manipulation.
50
51 For example, the FEM framework does not specify the number of spatial
52 dimensions.  Node locations are treated as just another kind of node data,
53 with no restrictions on the number of data items.
54 This allows the FEM framework to work with problems having any number 
55 of spatial dimensions.
56
57
58 \subsection{Terminology}
59 \label{sec:terminology}
60
61 A FEM program manipulates elements and nodes. An \term{element} is a portion of
62 the problem domain, also known as a cell, and is typically some simple shape 
63 like a triangle, square, or hexagon in 2D; or tetrahedron or rectangular solid in 3D.  
64 A \term{node} is a point in the domain, and is often the vertex of several elements.  
65 Together, the elements and nodes form a \term{mesh}, which is the 
66 central data structure in the FEM framework.
67
68 An element knows which nodes surround it via the element's
69 \term{connectivity table}, which lists the nodes adjacent to each element.
70
71 \begin{figure}[h]
72 \begin{center}
73 \includegraphics[width=4in]{fig/simple_mesh}
74 \end{center}
75 \caption{3-element, 5 node mesh.}
76 \label{fig:simplemesh}
77 \end{figure}
78
79 \begin{table}[h]
80 \begin{center}
81 \begin{tabular}{||l||l|l|l||}\hline
82 Element & \multicolumn{3}{c||}{Adjacent Nodes} \\\hline
83 e1 & n1 & n3 & n4 \\
84 e2 & n1 & n2 & n4 \\
85 e3 & n2 & n4 & n5 \\
86 \hline
87 \end{tabular}
88 \end{center}
89 \caption{Connectivity table for mesh in figure~\ref{fig:simplemesh}.}
90 \label{table:simplemesh}
91 \end{table}
92
93 A typical FEM program performs some element-by-element calculations which
94 update adjacent node values; then some node-by-node calculations.  For
95 example, a material dynamics program has the structure:
96
97 \begin{alltt}
98      time loop
99           element loop-- Element deformation applies forces to
100           surrounding nodes
101           node loop-- Forces and boundary conditions change node
102           positions
103      end time loop
104 \end{alltt}
105
106 We can parallelize such FEM programs by partitioning the serial mesh
107 elements into several smaller meshes, or \term{chunks}.  There is normally
108 at least one chunk per processor; and often even more.  During partitioning, 
109 we give nodes and elements new, \term{local} numbers within that chunk.
110 In the figure below, we have partitioned the mesh above into two chunks, A and B.
111
112 \begin{figure}[h]
113 \begin{center}
114 \includegraphics[width=4in]{fig/partitioned_mesh}
115 \end{center}
116 \caption{Partitioned mesh.}
117 \label{fig:partitionedmesh}
118 \end{figure}
119
120 \begin{table}[hh]
121 \begin{center}
122 \begin{tabular}{||l||l|l|l||}\hline
123 Element & \multicolumn{3}{c||}{Adjacent Nodes} \\\hline
124 e1 & n1 & n3 & n4 \\
125 e2 & n1 & n2 & n4 \\
126 \hline
127 \end{tabular}
128 \end{center}
129 \caption{Connectivity table for chunk A in figure~\ref{fig:partitionedmesh}.}
130 \label{table:chunkA}
131 \end{table}
132
133 \begin{table}[hh]
134 \begin{center}
135 \begin{tabular}{||l||l|l|l||}\hline
136 Element & \multicolumn{3}{c||}{Adjacent Nodes}\\\hline
137 e1 & n1 & n2 & n3 \\
138 \hline
139 \end{tabular}
140 \end{center}
141 \caption{Connectivity table for chunk B in figure~\ref{fig:partitionedmesh}.}
142 \label{table:chunkB}
143 \end{table}
144
145 Note that chunk A's node n2 and B's node n1 were actually the same node in
146 the original mesh-- partitioning split this single node into two shared
147 copies (one on each chunk).  However, since adding forces is associative, we
148 can handle shared nodes by computing the forces normally (ignoring the
149 existence of the other chunk), then adding both chunks' net force for the
150 shared node together.  This ``node update'' will give us the same resulting
151 force on each shared node as we would get without partitioning, thus the
152 same positions, thus the same final result.  
153
154 For example, under hydrostatic pressure, each chunk might compute a local
155 net force vector for its nodes as shown in Figure~\ref{fig:forcedecomp}
156 (a).  After adding forces across chunks, we have the consistent global forces
157 shown in Figure~\ref{fig:forcedecomp} (b).
158
159 \begin{figure}[h]
160 \begin{center}
161 \includegraphics[height=3in]{fig/forcedecomp}
162 \end{center}
163 \caption{A force calculation decomposed across chunks: (a) before update
164 (b) after updating forces across nodes.}
165 \label{fig:forcedecomp}
166 \end{figure}
167
168 Hence, each chunk's time loop has the structure:
169
170 \begin{alltt}
171      chunk time loop
172           element loop-- Element deformation applies forces to
173           surrounding nodes
174           <update forces on shared nodes>
175           node loop-- Forces and boundary conditions change node
176           positions
177      end time loop
178 \end{alltt}
179
180 This is exactly the form of the time loop for a \charmpp{} FEM framework
181 program.  The framework will accept a serial mesh, partition it, distribute
182 the chunks to each processor, then you run your time loop to perform
183 analysis and communication.
184
185
186 \subsection{Structure of a Classic FEM Framework Program}
187
188 A classic FEM framework program consists of two subroutines: 
189 \kw{init()} and \kw{driver()}.
190 \kw{init()} is called by the FEM framework
191 only on the first processor -- this routine typically does specialized I/O,
192 startup and shutdown tasks.  \kw{driver()} is called for every chunk on every
193 processor, and does the main work of the program.  In the language of the
194 TCHARM manual, \kw{init()} runs in the serial context, 
195 and \kw{driver()} runs in the parallel context.
196
197 \begin{alltt}
198      subroutine init
199           read the serial mesh and configuration data
200      end subroutine
201 /* after init, the FEM framework partitions the mesh */
202      subroutine driver
203           get local mesh chunk
204           time loop
205                FEM computations
206                communicate boundary conditions
207                more FEM computations
208           end time loop
209      end subroutine
210 \end{alltt}
211
212 In this mode, the FEM framework sets up a default writing
213 mesh during \kw{init()}, partitions the mesh after \kw{init()},
214 and sets up the partitioned mesh as the default reading mesh 
215 during \kw{driver()}. 
216
217
218 \subsection{Structure of an AMPI FEM Framework Program}
219
220 In addition to the classic init/driver structure above,
221 you can write an FEM framework program using the MPI style.
222 This is a more general, more flexible method of running
223 the program, but it is more complicated than the classic mode.
224 All FEM framework calls are available in either mode.
225
226 \begin{alltt}
227    main program
228       MPI_Init
229       FEM_Init(MPI_COMM_WORLD)
230       if (I am master processor)
231          read mesh
232       partition mesh
233       time loop
234           FEM computations
235           communicate boundary conditions
236           more FEM computations
237       end time loop
238    end main program
239 \end{alltt}
240
241 In this mode, the FEM framework does not set a default
242 reading or writing mesh, and does no partitioning;
243 so you must use the FEM\_Mesh routines to create and 
244 partition your mesh.
245 See the AMPI manual for details on how to declare
246 the main routine.
247
248 The driver() portion of a classic FEM program
249 strongly resembles an MPI mode main routine---in fact, a classic
250 FEM program can even make MPI calls from its \kw{driver()} 
251 routine, because the FEM framework is implemented directly on
252 top of MPI.  
253
254 There is even a special shell script for collecting
255 up the FEM framework source code to build a non-Charm, 
256 MPI-only version of the FEM framework.
257 To build FEM in this manner, you first build Charm++ normally,
258 then run a script to collect up the neccessary source files
259 (the FEM framework, a small number of Charm configuration 
260 and utility files, and the METIS library),
261 and finally build the library using the usual MPI compiler 
262 commands:
263 \begin{alltt}
264  > cd charm/
265  > ./src/libs/ck-libs/fem/make_fem_alone.sh
266  > cd fem_alone/
267  > mpicc -I. -DFEM_ALONE=1 -c *.c *.C 
268  > ar cr libfem_alone.a *.o
269 \end{alltt}
270 You will then have to build your application with the MPI
271 compilers, and manually point to this ``fem\_alone'' 
272 directory to find include files and the new FEM library.
273 A typical compiler invocation would be:
274 \begin{alltt}
275  > mpif90 -I$HOME/charm/fem_alone -L$HOME/charm/fem_alone foo.f90 -lfem_alone -o foo
276 \end{alltt}
277 This ``standalone'', non-Charm++ method of building the 
278 FEM framework prevents the use of load balancing or the other
279 features of Charm++, so we do not recommend it for normal use.
280
281
282 \subsection{Compilation and Execution}
283
284 A FEM framework program is a \charmpp\ program, so you must begin by
285 downloading the latest source version of \charmpp\ from
286 {\tt http://charm.cs.uiuc.edu/}.  Build the source with 
287 {\tt ./build FEM version} or {\tt cd} into the build directory, 
288 {\tt version/tmp}, and type {\tt make FEM}.
289 To compile a FEM program, pass the {\tt -language fem} (for C) or 
290 {\tt -language femf} (for Fortran) option to {\tt charmc}.
291 You can also build using the ``fem\_alone'' mode described
292 at the end of the section above.
293
294 In a charm installation, see charm/version/pgms/charm++/fem/
295 for several example and test programs.
296
297 At runtime, a Charm++/FEM framework program accepts the following
298 options, in addition to all the usual Charm++ options described in 
299 the Charm++ ``Installation and Usage Manual''.
300
301 \begin{itemize}
302 \item {\tt +vp} $v$  
303
304 Create $v$ mesh chunks, or ``virtual processors''.
305 By default, the number of mesh chunks is equal to the number of 
306 physical processors (set with {\tt +p} $p$).
307
308
309 \item {\tt -write}
310
311 Skip \kw{driver()}.
312 After running \kw{init()} normally, the framework partitions the mesh, 
313 writes the mesh partitions to files, and exits.  As usual, the
314 {\tt +vp} $v$ option controls the number of mesh partitions.
315
316 This option is only used in the classic mode---MPI-style programs
317 are not affected.
318
319
320 \item {\tt -read}
321
322 Skip \kw{init()}.
323 The framework reads the partitioned input mesh from files
324 and calls \kw{driver()}.  Together with {\tt -write}, this option
325 allows you to separate out the mesh preparation and partitioning 
326 phase from the actual parallel solution run.
327
328 This can be useful, for example, if \kw{init()} requires more memory 
329 to hold the unpartitioned mesh than is available on one processor of 
330 the parallel machine.  To avoid this limitation, you can run the program
331 with {\tt -write} on a machine with a lot of memory to prepare the input
332 files, then copy the files and run with {\tt -read} on a machine with 
333 a lot of processors.
334
335 {\tt -read} can also be useful during debugging or performance tuning, 
336 by skipping the (potentially slow) mesh preparation phase.
337 This option is only used in the classic mode---MPI-style programs
338 are not affected.
339
340
341 \item {\tt +tcharm\_trace fem}
342
343 Give a diagnostic printout on every call into the FEM framework.
344 This can be useful for locating a sudden crash, or understanding
345 how the program and framework interact.  Because printing the 
346 diagnostics can slow a program down, use this option with care.
347
348 \end{itemize}
349
350
351 \section{FEM Framework API Reference}
352
353 Some of the routines in the FEM framework have different requirements or meanings
354 depending on where they are called from.  When a routine is described
355 as being ``called from driver'', this means it is called in the parallel
356 context---from \kw{driver()} itself, any subroutine called by \kw{driver()},
357 or from whatever routine is run by the FEM-attached TCHARM threads.
358 When a routine is described as being ``called from init'', this means it is 
359 called in the serial context---from \kw{init()} itself, from any subroutine
360 called from \kw{init()}, from a routine called by \kw{FEM\_Update\_mesh},
361 or from whatever TCHARM code executes before the \kw{FEM\_Attach}.
362
363
364 \subsection{Utility}
365
366 \prototype{FEM\_Num\_partitions}
367 \function{int FEM\_Num\_partitions();}
368 \function{INTEGER FUNCTION :: FEM\_Num\_partitions()}
369
370      Return the number of mesh chunks in the current computation.  Can
371      only be called from the driver routine.
372
373 \prototype{FEM\_My\_partitions}
374 \function{int FEM\_My\_partition();}
375 \function{INTEGER FUNCTION :: FEM\_My\_partition()}
376
377      Return the number of the current chunk, from 0 to
378      \kw{num\_partitions}-1.  Can only be called from the driver routine.
379
380 \prototype{FEM\_Timer}
381 \function{double FEM\_Timer();}
382 \function{DOUBLE PRECISION FUNCTION :: FEM\_Timer()}
383
384      Return the current wall clock time, in seconds.  Resolution is
385      machine-dependent, but is at worst 10ms.
386
387 \prototype{FEM\_Print\_partition}
388 \function{void FEM\_Print\_partition();}
389 \function{SUBROUTINE FEM\_Print\_partition()}
390
391      Print a debugging representation of the current chunk's mesh.
392      Prints the entire connectivity array, and data associated with
393      each local node and element.
394
395 \prototype{FEM\_Print}
396 \function{void FEM\_Print(const char *str);}
397 \function{SUBROUTINE FEM\_Print(str)}
398 \args{CHARACTER*, INTENT(IN) :: str}
399
400      Print the given string, with "[<chunk number>]" printed 
401      before the text.  
402
403      This routine is no longer required: you can now use 
404      the usual printf, PRINT, or WRITE statements.
405
406
407 \input{mesh}
408
409 \input{idxl}
410
411 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
412 \section{Old Communication Routines}
413
414 (This section is for backward compatability only.  The IDXL routines
415 are the new, more flexible way to perform communication.)
416
417 The FEM framework handles the updating of the values of shared nodes-- that
418 is, it combines shared nodes' values across all processors.  The basic
419 mechanism to do this update is the ``field''-- numeric data items associated
420 with each node. We make no assumptions about the meaning of the node data,
421 allow various data types, and allow a mix of communicated and non-communicated 
422 data associated with each node.  The framework uses IDXL layouts to find the data items
423 associated with each node in memory.
424
425 Each field represents a (set of) data records stored in a contiguous array,
426 often indexed by node number.  You create a field once, with the IDXL layout 
427 routines or \kw{FEM\_Create\_field},
428 then pass the resulting field ID to \kw{FEM\_Update\_field} (which does the
429 shared node communication), \kw{FEM\_Reduce\_field} (which applies a
430 reduction over node values), or one of the other routines described below.
431
432
433 \prototype{FEM\_Update\_field}
434 \function{void FEM\_Update\_field(int Fid,void *nodes);}
435 \function{SUBROUTINE FEM\_Update\_field(Fid,nodes)}
436   \args{INTEGER, INTENT(IN)  :: Fid}
437   \args{varies, INTENT(INOUT) :: nodes}
438
439      Combine a field of all shared nodes with the other chunks.  Sums
440      the value of the given field across all chunks that share each
441      node.  For the example above, once each chunk has computed the net
442      force on each local node, this routine will sum the net force
443      across all shared nodes.
444
445      \kw{FEM\_Update\_field} can only be called from driver, and to be useful,
446      must be called from every chunk's driver routine.
447
448      After this routine returns, the given field of each shared node
449      will be the same across all processors that share the node.
450      
451      This routine is eqivalent to an \kw{IDXL\_Comm\_Sendsum} operation.
452
453 \prototype{FEM\_Read\_field}
454 \function{void FEM\_Read\_field(int Fid,void *nodes,char *fName);}
455 \function{SUBROUTINE FEM\_Read\_field(Fid,nodes,fName)}
456   \args{INTEGER, INTENT(IN)  :: Fid}
457   \args{varies, INTENT(OUT) :: nodes}
458   \args{CHARACTER*, INTENT(IN) :: fName}
459
460      Read a field out of the given serial input file.  The serial input
461      file is line-oriented ASCII-- each line begins with the global
462      node number (which must match the line order in the file),
463      followed by the data to be read into the node field.  The
464      remainder of each line is unread.  If called from Fortran, the
465      first line must be numbered 1; if called from C, the first line
466      must be numbered zero.  All fields are separated by white space
467      (any number of tabs or spaces).
468
469      For example, if we have called \kw{Create\_field} to describe 3 doubles,
470      the input file could begin with
471
472 \begin{alltt}
473           1    0.2    0.7    -0.3      First node
474           2    0.4    1.12   -17.26    another node
475           ...
476 \end{alltt}
477
478      \kw{FEM\_Read\_field} must be called from driver at any time, independent
479      of other chunks. 
480      
481      This routine has no IDXL equivalent.
482
483 \prototype{FEM\_Reduce\_field}
484 \function{void FEM\_Reduce\_field(int Fid,const void *nodes,void *out,int op);}
485 \function{SUBROUTINE FEM\_Reduce\_field(Fid,nodes,outVal,op)}
486   \args{INTEGER, INTENT(IN)  :: Fid,op}
487   \args{varies, INTENT(IN) :: nodes}
488   \args{varies, INTENT(OUT) :: outVal}
489
490 Combine one record per node of this field, according to op, across all chunks.
491 Shared nodes are not double-counted-- only one copy will contribute to the
492 reduction.  After \kw{Reduce\_field} returns, all chunks will have identical
493 values in \kw{outVal,} which must be \kw{vec\_len} copies of \kw{base\_type}.
494
495      May only be called from driver, and to complete, must be called
496      from every chunk's driver routine.
497
498      \kw{op} must be one of:
499
500 \begin{itemize}
501         \item \kw{FEM\_SUM}-- each element of \kw{outVal} will be the sum 
502 of the corresponding fields of all nodes
503         \item \kw{FEM\_MIN}-- each element of \kw{outVal} will be the 
504 smallest value among the corresponding field of all nodes
505         \item \kw{FEM\_MAX}-- each element of \kw{outVal} will be the largest 
506 value among the corresponding field of all nodes
507 \end{itemize}
508
509      This routine has no IDXL equivalent.
510
511
512 \prototype{FEM\_Reduce}
513 \function{void FEM\_Reduce(int Fid,const void *inVal,void *outVal,int op);}
514 \function{SUBROUTINE FEM\_Reduce(Fid,inVal,outVal,op)}
515   \args{INTEGER, INTENT(IN)  :: Fid,op}
516   \args{varies, INTENT(IN) :: inVal}
517   \args{varies, INTENT(OUT) :: outVal}
518
519      Combine one record of this field from each chunk, acoording to \kw{op}, 
520  across all chunks.
521 \kw{Fid} is only used for the \kw{base\_type} and \kw{vec\_len}-- offset and
522 \kw{dist} are not used.  After this call returns, all chunks will have
523 identical values in \kw{outVal}.  Op has the same values and meaning as
524 \kw{FEM\_Reduce\_field}.
525
526      May only be called from driver, and to complete, must be called
527      from every chunk's driver routine.
528
529 \begin{alltt}
530 ! C example
531    double inArr[3], outArr[3];
532    int fid=IDXL_Layout_create(FEM_DOUBLE,3);
533    FEM_Reduce(fid,inArr,outArr,FEM_SUM);
534
535 ! f90 example
536    DOUBLE PRECISION :: inArr(3), outArr(3)
537    INTEGER fid
538    fid=IDXL_Layout_create(FEM_DOUBLE,3)
539    CALL FEM_Reduce(fid,inArr,outArr,FEM_SUM)
540 \end{alltt}
541
542      This routine has no IDXL equivalent.
543
544
545 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
546 \subsection{Ghost Communication}
547
548 It is possible to get values for a chunk's ghost nodes and elements from the neighbors. To do this, use:
549
550 \prototype{FEM\_Update\_ghost\_field}
551 \function{void FEM\_Update\_ghost\_field(int Fid, int elTypeOrMinusOne, void *data);}
552 \function{SUBROUTINE FEM\_Update\_ghost\_field(Fid,elTypeOrZero,data)}
553   \args{INTEGER, INTENT(IN)  :: Fid,elTypeOrZero}
554   \args{varies, INTENT(INOUT) :: data}
555
556 This has the same requirements and call sequence as \kw{FEM\_Update\_field}, except it applies to ghosts. You specify which type of element to exchange using the elType parameter. Specify -1 (C version) or 0 (fortran version) to exchange node values.  
557
558
559 \subsection{Ghost List Exchange}
560
561 It is possible to exchange sparse lists of ghost elements between FEM chunks.
562
563 \prototype{FEM\_Exchange\_ghost\_lists}
564 \function{void FEM\_Exchange\_ghost\_lists(int elemType,int nIdx,const int *localIdx);}
565 \function{SUBROUTINE FEM\_Exchange\_ghost\_lists(elemType,nIdx,localIdx)}
566   \args{INTEGER, INTENT(IN)  :: elemType,nIdx}
567   \args{INTEGER, INTENT(IN) :: localIdx[nIdx]}
568
569 This routine sends the local element indices in localIdx to those neighboring chunks that connect to its ghost elements on the other side.  That is, if the element
570 \kw{localIdx[i]} has a ghost on some chunk \kw{c}, \kw{localIdx[i]} will be sent to 
571 and show up in the ghost list of chunk \kw{c}.
572
573 \prototype{FEM\_Get\_ghost\_list\_length}
574 \function{int FEM\_Get\_ghost\_list\_length(void);}
575 Returns the number of entries in my ghost list---the number of my ghosts that
576 other chunks passed to their call to \kw{FEM\_Exchange\_ghost\_lists}.
577
578 \prototype{FEM\_Get\_ghost\_list}
579 \function{void FEM\_Get\_ghost\_list(int *retLocalIdx);}
580 \function{SUBROUTINE FEM\_Get\_ghost\_list(retLocalIdx)}
581   \args{INTEGER, INTENT(OUT) :: retLocalIdx[FEM\_Get\_ghost\_list\_length()]}
582
583 These routines access the list of local elements sent by other chunks.  
584 The returned indices will all refer to ghost elements in my chunk.
585
586
587
588 % Permanently undocumented (i.e., officially denied) features:
589 %   FEM_Serial_Split(ndom) ! Partition into ndom domains
590 %   FEM_Serial_Begin(idom) ! Begin accessing the idom'th domain
591 %   
592 % ``There has to be a better way:''
593 %   FEM_Composite_elem, FEM_Combine_elem
594
595
596 \input{parfum}
597
598 \input{index}
599 \end{document}