c0cbcc840dda048ccc1172313f62eba946b8ac94
[charm.git] / doc / fem / manual.rst
1 =====================================
2 Finite Element Method (FEM) Framework
3 =====================================
4
5 .. contents::
6    :depth: 3
7
8 Introduction
9 ============
10
11 The Finite Element Method (FEM) approach is used in many engineering
12 applications with irregular domains, from elastic deformation problems
13 to crack propagation to fluid flow. Charm++ is a free message-passing
14 parallel runtime system for machines from clusters of workstations to
15 tightly-coupled SMPs. The Charm++ FEM framework allows you to write a
16 parallel FEM program, in C or Fortran 90, that closely resembles a
17 serial version but includes a few framework calls.
18
19 Using the FEM framework also allows you to take advantage of all the
20 features of Charm++, including run-time load balancing, performance
21 monitoring and visualization, and checkpoint/restart, with no additional
22 effort. The FEM framework also combines naturally with other Charm++
23 frameworks built on TCHARM.
24
25 The FEM framework has been undergoing a wave of recent improvements. A
26 choice to rename the new version ParFUM has been adopted. ParFUM is short
27 for Parallel Framework for Unstructured Meshes. Section
28 :numref:`sec:ParFUM` describes some of the new features included in
29 ParFUM that were not present in FEM.
30
31 Philosophy
32 ----------
33
34 The Charm++ FEM framework is designed to be flexible, in that it
35 provided a few very general operations, such as loading and partitioning
36 a “mesh.” In describing these operations, we draw on examples from
37 structural analysis, but in fact the same calls can be used for other
38 applications, including fluid dynamics or partial differential equations
39 solvers, or even general-purpose graph manipulation.
40
41 For example, the FEM framework does not specify the number of spatial
42 dimensions. Node locations are treated as just another kind of node
43 data, with no restrictions on the number of data items. This allows the
44 FEM framework to work with problems having any number of spatial
45 dimensions.
46
47 .. _sec:terminology:
48
49 Terminology
50 -----------
51
52 A FEM program manipulates elements and nodes. An **element** is a
53 portion of the problem domain, also known as a cell, and is typically
54 some simple shape like a triangle, square, or hexagon in 2D; or
55 tetrahedron or rectangular solid in 3D. A **node** is a point in the
56 domain, and is often the vertex of several elements. Together, the
57 elements and nodes form a **mesh**, which is the central data structure
58 in the FEM framework.
59
60 An element knows which nodes surround it via the element’s
61 **connectivity table**, which lists the nodes adjacent to each element.
62
63 .. figure:: fig/simple_mesh.png
64    :name: fig:simplemesh
65    :width: 4in
66
67    3-element, 5 node mesh.
68
69 .. table:: Connectivity table for mesh in figure :numref:`fig:simplemesh`.
70
71    ======= ==== ==== ====
72    Element Adjacent Nodes
73    ======= ==============
74    e1      n1   n3   n4
75    e2      n1   n2   n4
76    e3      n2   n4   n5
77    ======= ==== ==== ====
78
79 A typical FEM program performs some element-by-element calculations
80 which update adjacent node values; then some node-by-node calculations.
81 For example, a material dynamics program has the structure:
82
83 .. code-block:: none
84
85         time loop
86              element loop-- Element deformation applies forces to
87              surrounding nodes
88              node loop-- Forces and boundary conditions change node
89              positions
90         end time loop
91
92 We can parallelize such FEM programs by partitioning the serial mesh
93 elements into several smaller meshes, or **chunks**. There is normally
94 at least one chunk per processor; and often even more. During
95 partitioning, we give nodes and elements new, **local** numbers within
96 that chunk. In the figure below, we have partitioned the mesh above into
97 two chunks, A and B.
98
99 .. figure:: fig/partitioned_mesh.png
100    :name: fig:partitionedmesh
101    :width: 4in
102
103    Partitioned mesh.
104
105 .. table:: Connectivity table for chunk A in figure :numref:`fig:partitionedmesh`.
106
107    ======= ==== ==== ====
108    Element Adjacent Nodes
109    ======= ==============
110    e1      n1   n3   n4
111    e2      n1   n2   n4
112    ======= ==== ==== ====
113
114 .. table:: Connectivity table for chunk B in figure :numref:`fig:partitionedmesh`.
115
116    ======= ==== ==== ====
117    Element Adjacent Nodes
118    ======= ==============
119    e1      n1   n2   n3
120    ======= ==== ==== ====
121
122 Note that chunk A’s node n2 and B’s node n1 were actually the same node
123 in the original mesh- partitioning split this single node into two
124 shared copies (one on each chunk). However, since adding forces is
125 associative, we can handle shared nodes by computing the forces normally
126 (ignoring the existence of the other chunk), then adding both chunks’
127 net force for the shared node together. This “node update” will give us
128 the same resulting force on each shared node as we would get without
129 partitioning, thus the same positions, thus the same final result.
130
131 For example, under hydrostatic pressure, each chunk might compute a
132 local net force vector for its nodes as shown in
133 Figure :numref:`fig:forcedecomp` (a). After adding forces across
134 chunks, we have the consistent global forces shown in
135 Figure :numref:`fig:forcedecomp` (b).
136
137 .. figure:: fig/forcedecomp.png
138    :name: fig:forcedecomp
139    :height: 3in
140
141    A force calculation decomposed across chunks: (a) before update (b)
142    after updating forces across nodes.
143
144 Hence, each chunk’s time loop has the structure:
145
146 ::
147
148         chunk time loop
149              element loop-- Element deformation applies forces to
150              surrounding nodes
151              <update forces on shared nodes>
152              node loop-- Forces and boundary conditions change node
153              positions
154         end time loop
155
156 This is exactly the form of the time loop for a Charm++ FEM framework
157 program. The framework will accept a serial mesh, partition it,
158 distribute the chunks to each processor, then you run your time loop to
159 perform analysis and communication.
160
161 Structure of a Classic FEM Framework Program
162 --------------------------------------------
163
164 A classic FEM framework program consists of two subroutines: init() and
165 driver(). init() is called by the FEM framework only on the first
166 processor - this routine typically does specialized I/O, startup and
167 shutdown tasks. driver() is called for every chunk on every processor,
168 and does the main work of the program. In the language of the TCHARM
169 manual, init() runs in the serial context, and driver() runs in the
170 parallel context.
171
172 ::
173
174         subroutine init
175              read the serial mesh and configuration data
176         end subroutine
177    /* after init, the FEM framework partitions the mesh */
178         subroutine driver
179              get local mesh chunk
180              time loop
181                   FEM computations
182                   communicate boundary conditions
183                   more FEM computations
184              end time loop
185         end subroutine
186
187 In this mode, the FEM framework sets up a default writing mesh during
188 init(), partitions the mesh after init(), and sets up the partitioned
189 mesh as the default reading mesh during driver().
190
191 Structure of an AMPI FEM Framework Program
192 ------------------------------------------
193
194 In addition to the classic init/driver structure above, you can write an
195 FEM framework program using the MPI style. This is a more general, more
196 flexible method of running the program, but it is more complicated than
197 the classic mode. All FEM framework calls are available in either mode.
198
199 ::
200
201       main program
202          MPI_Init
203          FEM_Init(MPI_COMM_WORLD)
204          if (I am master processor)
205             read mesh
206          partition mesh
207          time loop
208              FEM computations
209              communicate boundary conditions
210              more FEM computations
211          end time loop
212       end main program
213
214 In this mode, the FEM framework does not set a default reading or
215 writing mesh, and does no partitioning; so you must use the FEM_Mesh
216 routines to create and partition your mesh. See the AMPI manual for
217 details on how to declare the main routine.
218
219 The driver() portion of a classic FEM program strongly resembles an MPI
220 mode main routine—in fact, a classic FEM program can even make MPI calls
221 from its driver() routine, because the FEM framework is implemented
222 directly on top of MPI.
223
224 There is even a special shell script for collecting up the FEM framework
225 source code to build a non-Charm, MPI-only version of the FEM framework.
226 To build FEM in this manner, you first build Charm++ normally, then run
227 a script to collect up the neccessary source files (the FEM framework, a
228 small number of Charm configuration and utility files, and the METIS
229 library), and finally build the library using the usual MPI compiler
230 commands:
231
232 .. code-block:: bash
233
234    $ cd charm/
235    $ ./src/libs/ck-libs/fem/make_fem_alone.sh
236    $ cd fem_alone/
237    $ mpicc -I. -DFEM_ALONE=1 -c *.c *.C
238    $ ar cr libfem_alone.a *.o
239
240 You will then have to build your application with the MPI compilers, and
241 manually point to this “fem_alone” directory to find include files and
242 the new FEM library. A typical compiler invocation would be:
243
244 .. code-block:: bash
245
246    $ mpif90 -I$HOME/charm/fem_alone -L$HOME/charm/fem_alone foo.f90 -lfem_alone -o foo
247
248 This “standalone”, non-Charm++ method of building the FEM framework
249 prevents the use of load balancing or the other features of Charm++, so
250 we do not recommend it for normal use.
251
252 Compilation and Execution
253 -------------------------
254
255 A FEM framework program is a Charm++ program, so you must begin by
256 downloading the latest source version of Charm++ from
257 ``http://charm.cs.uiuc.edu/``. Build the source with
258 ``./build FEM version`` or ``cd`` into the build directory,
259 ``version/tmp``, and type ``make FEM``. To compile a FEM program, pass
260 the ``-language fem`` (for C) or ``-language femf`` (for Fortran) option
261 to ``charmc``. You can also build using the “fem_alone” mode described
262 at the end of the section above.
263
264 In a charm installation, see charm/version/pgms/charm++/fem/ for several
265 example and test programs.
266
267 At runtime, a Charm++/FEM framework program accepts the following
268 options, in addition to all the usual Charm++ options described in the
269 Charm++ “Installation and Usage Manual”.
270
271 -  ``+vp`` :math:`v`
272
273    Create :math:`v` mesh chunks, or “virtual processors”. By default,
274    the number of mesh chunks is equal to the number of physical
275    processors (set with ``+p`` :math:`p`).
276
277 -  ``-write``
278
279    Skip driver(). After running init() normally, the framework
280    partitions the mesh, writes the mesh partitions to files, and exits.
281    As usual, the ``+vp`` :math:`v` option controls the number of mesh
282    partitions.
283
284    This option is only used in the classic mode—MPI-style programs are
285    not affected.
286
287 -  ``-read``
288
289    Skip init(). The framework reads the partitioned input mesh from
290    files and calls driver(). Together with ``-write``, this option
291    allows you to separate out the mesh preparation and partitioning
292    phase from the actual parallel solution run.
293
294    This can be useful, for example, if init() requires more memory to
295    hold the unpartitioned mesh than is available on one processor of the
296    parallel machine. To avoid this limitation, you can run the program
297    with ``-write`` on a machine with a lot of memory to prepare the
298    input files, then copy the files and run with ``-read`` on a machine
299    with a lot of processors.
300
301    ``-read`` can also be useful during debugging or performance tuning,
302    by skipping the (potentially slow) mesh preparation phase. This
303    option is only used in the classic mode—MPI-style programs are not
304    affected.
305
306 -  ``+tcharm_trace fem``
307
308    Give a diagnostic printout on every call into the FEM framework. This
309    can be useful for locating a sudden crash, or understanding how the
310    program and framework interact. Because printing the diagnostics can
311    slow a program down, use this option with care.
312
313 FEM Framework API Reference
314 ===========================
315
316 Some of the routines in the FEM framework have different requirements or
317 meanings depending on where they are called from. When a routine is
318 described as being “called from driver”, this means it is called in the
319 parallel context—from driver() itself, any subroutine called by
320 driver(), or from whatever routine is run by the FEM-attached TCHARM
321 threads. When a routine is described as being “called from init”, this
322 means it is called in the serial context—from init() itself, from any
323 subroutine called from init(), from a routine called by FEM_Update_mesh,
324 or from whatever TCHARM code executes before the FEM_Attach.
325
326 Utility
327 -------
328
329 ::
330
331   int FEM_Num_partitions();
332
333 .. code-block:: fortran
334
335   INTEGER FUNCTION :: FEM_Num_partitions()
336
337 Return the number of mesh chunks in the current computation. Can only be
338 called from the driver routine.
339
340 ::
341
342   int FEM_My_partition();
343
344 .. code-block:: fortran
345
346   INTEGER FUNCTION :: FEM_My_partition()
347
348 Return the number of the current chunk, from 0 to num_partitions-1. Can
349 only be called from the driver routine.
350
351 ::
352
353   double FEM_Timer();
354
355 .. code-block:: fortran
356
357   DOUBLE PRECISION FUNCTION :: FEM_Timer()
358
359 Return the current wall clock time, in seconds. Resolution is
360 machine-dependent, but is at worst 10ms.
361
362 ::
363
364   void FEM_Print_partition();
365
366 .. code-block:: fortran
367
368   SUBROUTINE FEM_Print_partition()
369
370 Print a debugging representation of the current chunk’s mesh. Prints the
371 entire connectivity array, and data associated with each local node and
372 element.
373
374 ::
375
376   void FEM_Print(const char *str);
377
378 .. code-block:: fortran
379
380   SUBROUTINE FEM_Print(str)
381   CHARACTER*, INTENT(IN) ::  str
382
383 Print the given string, with "[<chunk number>]" printed before the
384 text.
385
386 This routine is no longer required: you can now use the usual printf,
387 PRINT, or WRITE statements.
388
389 .. _sec:entities:
390
391 Mesh Nodes and Elements
392 =======================
393
394 These routines describe and retrieve the finite element mesh for this
395 computation. A **mesh**, from the framework’s perspective, is a list of
396 elements, nodes, and other data that describes the computational domain.
397 The FEM framework provides extensive support for creating, manipulating,
398 and partitioning meshes.
399
400 A **serial mesh** consists of a single large piece. It’s usually easiest
401 to read and write serial meshes to existing, non-parallel file formats,
402 and it can be easier to manipulate serial meshes. By contrast, a
403 **parallel mesh** consists of several pieces, called **chunks** or
404 partitions. Different processors can work on different pieces of a
405 parallel mesh, so most of the computation is done using parallel meshes.
406 A simple program might create or read in a single serial mesh in init,
407 get a local chunk of the partitioned [1]_ mesh in driver, and work on
408 that chunk for the rest of the program. A more complex program might set
409 an initial mesh in init; then get, work on, reassemble and repartition
410 the mesh several times in driver via FEM_Update_mesh.
411
412 Mesh Entity Types
413 -----------------
414
415 A mesh consists of **entities**, such as nodes and elements. Entities
416 always have a **local number**, which is just the entities’ current
417 index in its array. Entites may also have a **global number**, which is
418 the entity’s index in the unpartitioned serial mesh. Entities have data
419 values called **attributes**. For example, the location of each node
420 might be called the “location” attribute of the “node” entity type.
421 Attributes are always stored in regular arrays indexed by the entity’s
422 local number. This table lists the different attributes that can be read
423 or written for each type of entity.
424
425 A **shared entity** is a boundary entity that two or more chunks can
426 both update—currently, only nodes can be shared. Shared nodes are mixed
427 in with regular nodes, and the framework currently provides no way to
428 identify which nodes are shared.
429
430 A **ghost entity** is a boundary entity that is asymmetrically
431 shared—one side provides values for the ghost from one of its real
432 entities, and the other sides accept read-only copies of these values.
433 Ghosts are described in more detail in Section :numref:`sec:ghost`, and
434 can be accessed by adding the constant FEM_GHOST to the corresponding
435 real entity’s type.
436
437 The different kinds of entities are described in the following sections.
438
439 =============================== =========================================
440 Real Entity                     Ghost Entity
441 =============================== =========================================
442 FEM_NODE                        FEM_GHOST+FEM_NODE
443 FEM_ELEM+\ :math:`elType`       FEM_GHOST+FEM_ELEM+\ :math:`elType`
444 FEM_SPARSE+\ :math:`sparseType` FEM_GHOST+FEM_SPARSE+\ :math:`sparseType`
445 =============================== =========================================
446
447 Nodes
448 ~~~~~
449
450 FEM_NODE is the entity code for nodes, the simplest kind of entity. A
451 node is a single point in the domain, and elements are defined by their
452 nodes. Nodes can have the following attributes:
453
454 -  FEM_DATA+\ :math:`tag` Uninterpreted user data, which might include
455    material properties, boundary conditions, flags, etc. User data can
456    have any data type and width. :math:`tag` can be any number from 0 to
457    one billion—it allows you to register several data fields with a
458    single entity.
459
460 -  FEM_GLOBALNO Global node numbers. Always a 1-wide index type.
461
462 -  FEM_SYMMETRIES Symmetries that apply to this node. Always a 1-wide
463    FEM_BYTE.
464
465 -  FEM_NODE_PRIMARY Marker indicating that this chunk is responsible for
466    this node. Every node is primary in exactly one chunk. This attribute
467    is always a 1-wide FEM_BYTE containing 0 or 1.
468
469 Elements
470 ~~~~~~~~
471
472 FEM_ELEM+\ :math:`elType` is the entity code for one kind of element.
473 :math:`elType` is a small, user-defined value that uniquely identifies
474 this element type. Like nodes, elements can have the attributes
475 FEM_DATA+\ :math:`tag`, FEM_GLOBALNO, or FEM_SYMMETRIES; but every
476 element type must have this attribute:
477
478 -  FEM_CONN Lists the numbers of the nodes around this element. See the
479    description in the ghost section for special ghost connectivity.
480    Always an index type-FEM_INDEX_0 for C-style 0-based node indexing,
481    or FEM_INDEX_1 for Fortran-style 1-based node indexing.
482
483 Sparse Elements
484 ~~~~~~~~~~~~~~~
485
486 FEM_SPARSE+\ :math:`sparseType` is the entity code for one kind of
487 sparse element. Again, :math:`sparseType` is a small, user-defined
488 unique value. The only difference between ordinary elements and sparse
489 elements regards partitioning. Ignoring ghosts, ordinary elements are
490 never duplicated—each element is sent to its own chunk. Sparse elements
491 may be duplicated, and are always dependent on some other entity for
492 their partitioning. Sparse elements have all the attributes of ordinary
493 elements: FEM_DATA+\ :math:`tag`, FEM_GLOBALNO, FEM_SYMMETRIES, and
494 FEM_CONN, as well as the special attribute FEM_SPARSE_ELEM.
495
496 Without the FEM_SPARSE_ELEM attribute, a sparse element will be copied
497 to every chunk that contains all the sparse element’s nodes. This is
498 useful for things like node-associated boundary conditions, where the
499 sparse element connectivity might list the nodes with boundary
500 conditions, and the sparse element data might list the boundary
501 condition values.
502
503 The FEM_SPARSE_ELEM attribute lists the ordinary element each sparse
504 element should be partitioned with. This attribute consists of pairs
505 (:math:`elType`,\ :math:`elNum`), indicating that this sparse element
506 should be sent to wherever the :math:`elNum`\ ’th
507 FEM_ELEM+\ :math:`elType` is partitioned.
508
509 -  FEM_SPARSE_ELEM Lists the element we should be partitioned with. The
510    width of this attribute is always 2, and the data type must be an
511    index type-FEM_INDEX_0 or FEM_INDEX_1.
512
513 Mesh Entity Manipulation
514 ------------------------
515
516 ::
517
518   int FEM_Mesh_default_read(void);
519
520 .. code-block:: fortran
521
522   INTEGER function :: FEM_Mesh_default_read()
523
524 Return the default reading mesh. This routine is valid:
525
526 -  From driver(), to return the partitioned mesh.
527
528 -  During your FEM_Update_mesh routine, to return the assembled mesh.
529
530 -  Anytime after a call to FEM_Mesh_set_default_read.
531
532 ::
533
534   int FEM_Mesh_default_write(void);
535
536 .. code-block:: fortran
537
538   INTEGER function :: FEM_Mesh_default_write()
539
540 Return the default writing mesh. This routine is valid:
541
542 -  From init(), to change the new serial mesh.
543
544 -  From driver(), to change the new partitioned mesh.
545
546 -  During your FEM_Update_mesh routine, to change the new serial mesh.
547
548 -  Anytime after a call to FEM_Mesh_set_default_write.
549
550 ::
551
552   int FEM_Mesh_get_length(int mesh,int entity);
553
554 .. code-block:: fortran
555
556   INTEGER function :: FEM_Mesh_get_length(mesh,entity)
557
558 Return the number of entitys that exist in this mesh.
559
560 This call can be used with any entity. For example, to get the number of
561 nodes,
562
563 ::
564
565          nNodes=FEM_Mesh_get_length(mesh,FEM_NODE)
566
567 To get the number of ghost nodes,
568
569 ::
570
571          nGhostNodes=FEM_Mesh_get_length(mesh,FEM_GHOST+FEM_NODE)
572
573 To get the number of real elements of type 2,
574
575 ::
576
577         nElem=FEM_Mesh_get_length(mesh,FEM_ELEM+2)
578
579
580 ::
581
582   void FEM_Mesh_data(int mesh,int entity,int attr, void *data, int
583   first, int length, int datatype,int width);
584
585 .. code-block:: fortran
586
587   SUBROUTINE FEM_Mesh_data(mesh,entity,attr,data,first,length,datatype,width)
588   INTEGER, INTENT(IN) :: mesh,entity,attr,first,length,datatype,width
589   datatype, intent(inout) :: data(width,length)
590
591 This is the one routine for getting or setting entity’s attributes on
592 the mesh.
593
594 -  mesh A FEM mesh object. Depending on whether this is a reading or
595    writing mesh, this routine reads from or writes to the data array you
596    pass in.
597
598 -  entity A FEM entity code, for example FEM_NODE or
599    FEM_GHOST+FEM_ELEM+1.
600
601 -  attr A FEM attribute code, for example FEM_DATA+\ :math:`tag` or
602    FEM_CONN.
603
604 -  data The user data to get or set. Each row of this array consists of
605    width values, and contains the data values of the attribute for the
606    corresponding entity. This data must be formatted as one of:
607
608    ::
609
610             datatype :: data(width,length)
611             datatype :: data(width*length)
612
613 -  first The first entity to affect. In C, this is normally 0; in
614    Fortran, this is normally 1.
615
616 -  length The number of entities to affect. The entities affected are
617    thus those numbered from first to first+length-1. For now, length
618    must be either 1, to touch a single entity; or else the total number
619    of entities-that is, FEM_Mesh_get_length(mesh,entity).
620
621 -  datatype The data type stored in this attribute. This is one of the
622    standard FEM data types FEM_BYTE, FEM_INT, FEM_FLOAT, or FEM_DOUBLE;
623    or else the C-style 0-based index type FEM_INDEX_0 or the
624    Fortran-style 1-based index type FEM_INDEX_1. Alternatively, the
625    equivalent types IDXL_BYTE, IDXL_INT, IDXL_FLOAT, IDXL_DOUBLE,
626    IDXL_INDEX_0, or IDXL_INDEX_1 may be used.
627
628 -  width The number of data items per entity.
629
630 For example, to set the element connectivity, which is stored as 3
631 integer node indices in nodes, you would:
632
633 ::
634
635   /* C version */
636   int *nodes=new int[3*nElems];
637   ... fill out nodes ...
638   FEM_Mesh_data(mesh,FEM_ELEM+1,FEM_CONN, nodes, 0,nElems, FEM_INDEX_0, 3);
639   ... continue to use or delete nodes ...
640
641 .. code-block:: fortran
642
643   ! F90 version
644   ALLOCATE(nodes(3,nElems))
645   ... fill out nodes ...
646   CALL FEM_Mesh_data(mesh,FEM_ELEM+1,FEM_CONN, nodes, 1,nElems, FEM_INDEX_1, 3)
647   ... continue to use or delete nodes ...
648
649 To add a new node property with 2 double-precision numbers from an array
650 mat (containing, for example, material properties), you would first pick
651 an unused user data "tag", for example 13, and:
652
653 ::
654
655   /* C version */
656   double *mat=new double[2*nNodes];
657   ...
658   FEM_Mesh_data(mesh,FEM_NODE, FEM_DATA+13, mat, 0,nNodes, FEM_DOUBLE, 2);
659
660 .. code-block:: fortran
661
662   ! F90 version
663   ALLOCATE(mat(2,nNodes))
664   CALL FEM_Mesh_data(mesh,FEM_NODE,FEM_DATA+13, mat, 1,nNodes, FEM_DOUBLE, 2)
665
666 Entity Inquiry
667 --------------
668
669 ::
670
671   int FEM_Mesh_get_width(int mesh,int entity,int attr);
672
673 .. code-block:: fortran
674
675   INTEGER function :: FEM_Mesh_get_width(mesh,entity,attr)
676   INTEGER, INTENT(IN) :: mesh,entity,attr
677
678 Return the width of the attribute attr of entity of mesh. This is the
679 value previously passed as “width” to FEM_Mesh_data.
680
681 ::
682
683   int FEM_Mesh_get_datatype(int mesh,int entity,int attr);
684
685 .. code-block:: fortran
686
687   INTEGER function :: FEM_Mesh_get_datatype(mesh,entity,attr)
688   INTEGER, INTENT(IN) :: mesh,entity,attr
689
690 Return the FEM data type of the attribute attr of entity of mesh. This
691 is the value previously passed as “datatype” to FEM_Mesh_data.
692
693 ::
694
695   int FEM_Mesh_get_entities(int mesh,int *entities);
696
697 .. code-block:: fortran
698
699   INTEGER function :: FEM_Mesh_get_entities(mesh,entities)
700   INTEGER, INTENT(IN) :: mesh
701   INTEGER, INTENT(OUT) :: entities(:)
702
703
704 Extract an array of the different entities present in this mesh.
705 Returns the number of entity types present. The entities array must be
706 big enough to hold all the different entities in the mesh.
707
708 For example, a simple mesh might have two entity types: FEM_NODE and
709 FEM_ELEM+1.
710
711
712 ::
713
714   int FEM_Mesh_get_attributes(int mesh,int entity,int *attributes);
715
716 .. code-block:: fortran
717
718   INTEGER function :: FEM_Mesh_get_attributes(mesh,entity,attributes)
719   INTEGER, INTENT(IN) :: mesh, entity
720   INTEGER, INTENT(OUT) :: attributes(:)
721
722 Extract an array of the different attributes of this entity. Returns
723 the number of attribute types present. The attributes array must be
724 big enough to hold all the attributes.
725
726 For example, a simple element might have three attributes: FEM_CONN for
727 node connectivity, FEM_GLOBALNO for global element numbers, and
728 FEM_DATA+7 for a material type.
729
730 ::
731
732   const char *FEM_Get_entity_name(int entity,char *storage);
733   const char *FEM_Get_attr_name(int attr,char *storage);
734   const char *FEM_Get_datatype_name(int datatype,char *storage);
735
736 Return a human-readable name for this FEM entity, attribute, or
737 datatype. The storage array must point to a buffer of at least 100
738 characters; this array might be used as temporary space to store the
739 returned string.
740
741 These routines are only available in C.
742
743 Advanced Entity Manipulation
744 ----------------------------
745 ::
746
747   void FEM_Mesh_data_offset(int mesh,int entity,int attr, void *data,
748   int first, int length, int datatype,int width, int offsetBytes,int
749   distanceBytes,int skewBytes);
750
751 .. code-block:: fortran
752
753   SUBROUTINE FEM_Mesh_data_offset(mesh,entity,attr,data,first,length,datatype,width,
754   offsetBytes,distanceBytes,skewBytes)
755   INTEGER, INTENT(IN) :: mesh,entity,attr,first,length,datatype,width
756   INTEGER, INTENT(IN) :: offsetBytes,distanceBytes,skewBytes
757   datatype, intent(inout) :: data(width,length)
758
759 This routine is a more complicated version of FEM_Mesh_data. It allows
760 you to get or set a mesh field directly from a user-defined structure.
761 See the documentation of IDXL_Layout_offset in
762 Section :numref:`sec:IDXLLayoutoffset` for details on how to set
763 offsetBytes, distanceBytes, and skewBytes.
764
765 ::
766
767   void FEM_Mesh_data_layout(int mesh,int entity,int attr, void *data,
768   int firstItem, int length, IDXL_Layout_t layout);
769
770 .. code-block:: fortran
771
772   SUBROUTINE FEM_Mesh_data_layout(mesh,entity,attr,data,first,length,layout)
773   INTEGER, INTENT(IN) :: mesh,entity,attr,first,length,layout
774   INTEGER, INTENT(IN) :: layout
775
776 This routine is a more complicated version of FEM_Mesh_data. Like
777 FEM_Mesh_data_offset, it allows you to get or set a mesh field
778 directly from a user-defined structure; but this routine expects the
779 structure to be described by an IDXL_Layout object.
780
781 .. _sec:mesh:
782
783 Meshes
784 ======
785
786 A "mesh" is a collection of nodes and elements knit together in memory,
787 as described in Section :numref:`sec:terminology`. Meshes are always
788 referred to by an integer that serves as a handle to the local mesh.
789
790 This section describes routines to manipulate entire meshes at once:
791 this includes calls to create and delete meshes, read and write meshes,
792 partition and reassemble meshes, and send meshes between processors.
793
794 Only a few of the mesh routines are collective; most of them only
795 describe local data and hence operate independently on each chunk.
796
797 Mesh Routines
798 -------------
799
800 ::
801
802   int FEM_Mesh_allocate(void);
803
804 .. code-block:: fortran
805
806   INTEGER FUNCTION :: FEM_Mesh_allocate()
807
808 Create a new local mesh object. The mesh is initially empty, but it is a
809 setting mesh, so call FEM_Mesh_data to fill the mesh with data.
810
811 ::
812
813   int FEM_Mesh_deallocate(int mesh);
814
815 .. code-block:: fortran
816
817   SUBROUTINE FEM_Mesh_deallocate(mesh)
818   INTEGER, INTENT(IN) :: mesh
819
820 Destroy this local mesh object, and its associated data.
821
822 ::
823
824   int FEM_Mesh_copy(int mesh);
825
826 .. code-block:: fortran
827
828   INTEGER FUNCTION FEM_Mesh_copy(mesh)
829   INTEGER, INTENT(IN) :: mesh
830
831 Create a new mesh object with a separate copy of the data stored in
832 this old mesh object.
833
834 ::
835
836   void FEM_Mesh_write(int mesh,const char *prefix,int partNo,int
837   nParts);
838
839 .. code-block:: fortran
840
841   SUBROUTINE FEM_Mesh_write(mesh,prefix,partNo,nParts)
842   INTEGER, INTENT(IN) :: mesh
843   INTEGER, INTENT(IN) :: partNo,nParts
844   character (LEN=*), INTENT(IN) :: prefix
845
846 Write this mesh to the file “prefix_vppartNo_nParts.dat”.
847
848 By convention, partNo begins at 0; but no index conversion is performed
849 so you can assign any meaning to partNo and nParts. In particular, this
850 routine is not collective-you can read any mesh from any processor. For
851 example, if prefix is “foo/bar”, the data for the first of 7 chunks
852 would be stored in “foo/bar_vp0_7.dat” and could be read using
853 FEM_Mesh_read(’foo/bar’,0,7).
854
855 Meshes are stored in a machine-portable format internal to FEM. The
856 format is currently ASCII based, but it is subject to change. We
857 strongly recommend using the FEM routines to read and write these files
858 rather than trying to prepare or parse them yourself.
859
860 ::
861
862   int FEM_Mesh_read(const char *prefix,int partNo,int nParts);
863
864 .. code-block:: fortran
865
866   INTEGER FUNCTION :: FEM_Mesh_read(prefix,partNo,nParts)
867   INTEGER, INTENT(IN) :: partNo,nParts
868   character (LEN=*), INTENT(IN) :: prefix
869
870 Read a new mesh from the file “prefix_vppartNo_nParts.dat”. The new
871 mesh begins in getting mode, so you can read the data out of the mesh
872 using calls to FEM_Mesh_data.
873
874 ::
875
876   int FEM_Mesh_broadcast(int mesh,int fromRank,FEM_Comm_t comm_context);
877
878 .. code-block:: fortran
879
880   INTEGER FUNCTION :: FEM_Mesh_broadcast(mesh,fromRank,comm_context)
881   INTEGER, INTENT(IN) :: mesh,fromRank,comm_context
882
883 Take the mesh mesh on processor fromRank (normally 0), partition the
884 mesh into one piece per processor (in the MPI communicator
885 comm_context, and return each processor its own piece of the
886 partitioned mesh. This call is collective, but only processor fromRank
887 needs to pass in a mesh; the mesh value is ignored on other
888 processors.
889
890 For example, if rank 0 has a mesh named “src”, we can partition src for
891 all the processors by executing:
892
893 ::
894
895      m=FEM_Mesh_broadcast(src,0,MPI_COMM_WORLD);
896
897 The new, partitioned mesh is in getting mode, so you can read the
898 partitioned data using calls to FEM_Mesh_data. This call does not affect
899 mesh in any way.
900
901 ::
902
903   int FEM_Mesh_reduce(int mesh,int toRank,FEM_Comm_t comm_context);
904
905 .. code-block:: fortran
906
907   INTEGER FUNCTION :: FEM_Mesh_reduce(mesh,toRank,comm_context)
908   INTEGER, INTENT(IN) :: mesh,toRank,comm_context
909
910 This call is the reverse operation of FEM_Mesh_broadcast: each
911 processor passes in a mesh in mesh, the mesh is assembled, and the
912 function returns the assembled mesh to processor toRank. This call is
913 collective, but only processor toRank is returned a mesh; all other
914 processors are returned the non-mesh value 0.
915
916 The new, reassembled mesh is in getting mode. This call does not affect
917 mesh.
918
919 Mesh Utility
920 ------------
921
922 ::
923
924   int FEM_Mesh_is_get(int mesh)
925
926 .. code-block:: fortran
927
928   INTEGER FUNCTION ::  FEM_Mesh_is_get(mesh)
929   INTEGER, INTENT(IN) :: mesh
930
931 Return true if this mesh is in getting mode. A getting mesh returns
932 values to FEM_Mesh_data.
933
934 ::
935
936   int FEM_Mesh_is_set(int mesh)
937
938 .. code-block:: fortran
939
940   INTEGER FUNCTION :: FEM_Mesh_is_set(mesh)
941   INTEGER, INTENT(IN) :: mesh
942
943 Return true if this mesh is in setting mode. A setting mesh extracts
944 values from FEM_Mesh_data.
945
946 ::
947
948   void FEM_Mesh_become_get(int mesh)
949
950 .. code-block:: fortran
951
952   SUBROUTINE :: FEM_Mesh_become_get(mesh)
953   INTEGER, INTENT(IN) :: mesh
954
955 Put this mesh in getting mode, so you can read back its values.
956
957 ::
958
959   void FEM_Mesh_become_set(int mesh)
960
961 .. code-block:: fortran
962
963   SUBROUTINE :: FEM_Mesh_become_set(mesh)
964   INTEGER, INTENT(IN) :: mesh
965
966 Put this mesh in setting mode, so you can set its values.
967
968 ::
969
970   void FEM_Mesh_print(int mesh);
971
972 .. code-block:: fortran
973
974   SUBROUTINE FEM_Mesh_print(mesh)
975   INTEGER, INTENT(IN) :: mesh
976
977 Print out a text description of the nodes and elements of this mesh.
978
979 Advanced Mesh Manipulation
980 --------------------------
981
982 ::
983
984   typedef void (*FEM_Userdata_fn)(pup_er p,void *data);
985   void FEM_Mesh_pup(int mesh,int pupTag, FEM_Userdata_fn fn, void *data);
986
987 .. code-block:: fortran
988
989   SUBROUTINE myPupFn(p,data)
990   INTEGER, INTENT(IN) :: p
991   TYPE(myType) :: data
992
993   SUBROUTINE FEM_Mesh_pup(mesh,pupTag,myPupFn,data)
994   INTEGER, INTENT(IN) :: mesh,pupTag
995   SUBROUTINE :: myPupFn
996   TYPE(myType) :: data
997
998
999 Store data with this mesh. data is a struct or TYPE with a pup
1000 function myPupFn—see the TCharm manual for details on writing a pup
1001 function. pupTag is an integer used to distinguish different pieces of
1002 data associated with this mesh.
1003
1004 When called on a setting mesh, this routine stores data; when called on
1005 a getting mesh, this routine reads out data.
1006
1007 data will be associated with the mesh itself, not any entity in the
1008 mesh. This makes it useful for storing shared data, often simulation
1009 constants such as the timestep or material properties. data is made a
1010 part of the mesh, and it will be read and written, sent and received,
1011 partitioned and assembled with the mesh.
1012
1013 ::
1014
1015   void FEM_Mesh_send(int mesh,int toRank,int tag,FEM_Comm_t
1016   comm_context);
1017
1018 .. code-block:: fortran
1019
1020   SUBROUTINE FEM_Mesh_send(mesh,toRank,tag,comm)
1021   INTEGER, INTENT(IN) :: mesh,toRank,tag,comm
1022
1023 Send the mesh mesh to the processor toRank, using the MPI tag tag and
1024 communicator comm_context. Tags are normally only needed if you plan
1025 to mix direct MPI calls with your FEM calls.
1026
1027 This call does not affect mesh.
1028
1029 ::
1030
1031   int FEM_Mesh_recv(int fromRank,int tag,FEM_Comm_t comm_context);
1032
1033 .. code-block:: fortran
1034
1035   INTEGER FUNCTION FEM_Mesh_recv(fromRank,tag,comm)
1036   INTEGER, INTENT(IN) :: fromRank,tag,comm
1037
1038 Receive a new mesh from the processor fromRank, using the MPI tag tag
1039 and communicator comm_context. You can also use the special values
1040 MPI_ANY_SOURCE as fromRank to receive a mesh from any processor, or
1041 use MPI_ANY_TAG for tag to match any tag.
1042
1043 The new mesh is returned in getting mode.
1044
1045 ::
1046
1047   void FEM_Mesh_partition(int mesh,int nParts,int *destMeshes);
1048
1049 .. code-block:: fortran
1050
1051   SUBROUTINE FEM_Mesh_partition(mesh,nParts,destMeshes)
1052   INTEGER, INTENT(IN) :: mesh,nParts
1053   INTEGER, INTENT(OUT) :: destMeshes(nParts)
1054
1055 Divide mesh into nParts pieces, and store the pieces into the array
1056 destMeshes.
1057
1058 The partitioned mesh is returned in getting mode. This is a local call;
1059 FEM_Mesh_broadcast is the collective version. This call does not affect
1060 the source mesh mesh.
1061
1062 ::
1063
1064   int FEM_Mesh_assemble(int nParts,const int *srcMeshes);
1065
1066 .. code-block:: fortran
1067
1068   INTEGER FUNCTION FEM_Mesh_assemble(nParts,srcMeshes)
1069   INTEGER, INTENT(IN) :: nParts, srcMeshes(nParts)
1070
1071 Assemble the nParts meshes listed in srcMeshes into a single mesh.
1072 Corresponding mesh pieces are matched using the attribute
1073 FEM_GLOBALNO. Specifically, if the value of the integer index
1074 attribute FEM_GLOBALNO for an entity is :math:`i`, the entity will be
1075 given the number :math:`i` in the reassembled mesh. If you do not set
1076 FEM_GLOBALNO, the different pieces of the mesh will remain
1077 separate—even “matching” nodes will not be merged.
1078
1079 The assembled mesh is returned in getting mode. This is a local call;
1080 FEM_Mesh_reduce is the collective version. This call does not affect the
1081 source meshes.
1082
1083 ::
1084
1085   void FEM_Mesh_copy_globalno(int src_mesh,int dest_mesh);
1086
1087 .. code-block:: fortran
1088
1089   SUBROUTINE FEM_Mesh_copy_globalno(src_mesh,dest_mesh)
1090   INTEGER, INTENT(IN) :: src_mesh,dest_mesh
1091
1092 Copy the FEM_GLOBALNO attribute for all the entity types in src_mesh
1093 into all the matching types in dest_mesh, where the matching types
1094 exist. This call is often used before an FEM_Mesh_assemble or
1095 FEM_Mesh_reduce to synchronize global numbers before reassembly.
1096
1097 .. _sec:ghost:
1098
1099 Mesh Ghosts
1100 ===========
1101
1102 A **ghost entity** is a local, read-only copy of a real entity on
1103 another chunk. Ghosts are typically added to the boundary of a chunk to
1104 allow the real (non-ghost) elements at the boundary to access values
1105 across the processor boundary. This makes a chunk “feel” as if it was
1106 part of a complete unpartitioned mesh; and can be useful with
1107 cell-centered methods, and in mesh modification.
1108
1109 .. figure:: fig/ghost_pre.png
1110    :name: fig:ghostpre
1111    :width: 1.5in
1112
1113    A small mesh partitioned into two pieces.
1114
1115 .. figure:: fig/ghost_edge.png
1116    :name: fig:ghostedge
1117    :width: 1.5in
1118
1119    The same mesh with one layer of edge-adjacent ghosts.
1120
1121 .. figure:: fig/ghost_node.png
1122    :name: fig:ghostnode
1123    :width: 1.5in
1124
1125    The same mesh with one layer of node-adjacent ghosts.
1126
1127 In Figure :numref:`fig:ghostpre`, we begin with a small
1128 mesh partitioned into pieces on the left and right. In
1129 Figure :numref:`fig:ghostedge`, we have added ghost
1130 elements (dark hashing) that share an edge with adjacent real elements
1131 (light hatching). In Figure :numref:`fig:ghostnode`, we add ghost
1132 elements that share at least one node with adjacent real elements.
1133
1134 .. _sec:ghostnum:
1135
1136 Ghost Numbering
1137 ---------------
1138
1139 Ghosts and real entities are stored by the framework in separate
1140 lists—to access the ghost entity type, add FEM_GHOST to the real
1141 entity’s type. For example, FEM_GHOST+FEM_ELEM+1 lists the ghost
1142 elements for elType 1. To get the number of ghost nodes, you would call
1143 FEM_Mesh_get_length(mesh,FEM_GHOST+FEM_NODE).
1144
1145 .. figure:: fig/conn_indexing.png
1146    :name: fig:connindexing
1147    :width: 4in
1148
1149    Node indices used in the element connectivity array. There are
1150    :math:`n` real nodes and :math:`m` ghosts.
1151
1152 For real elements, the element connectivity always consists of real
1153 nodes. But for ghost elements, the adjacent nodes may be missing, or may
1154 themselves be ghosts. Thus ghost element connectivity lists may include
1155 the invalid value -1 (in C) or 0 (in Fortran) to indicate that the
1156 corresponding node is not present; or may include values less than this,
1157 which indicate the corresponding node is a ghost. In C, ghost node
1158 :math:`i` is indicated by the value :math:`-2-i`, while in Fortran,
1159 ghost node :math:`i` is indicated by the value :math:`-i`. This node
1160 indexing system is illustrated in Figure :numref:`fig:connindexing`,
1161 This indexing system is bizarre, but it allows us to keep the real and
1162 ghost nodes clearly separate, while still allowing real and ghost nodes
1163 to be added in increasing order at both ends.
1164
1165 Since the C tests are complicated, in C we recommend using these macros:
1166
1167 -  FEM_Is_ghost_index(i) returns true if :math:`i` represents a ghost
1168    node. In Fortran, use the test :math:`i` .lt. :math:`0`
1169
1170 -  FEM_From_ghost_index(i) returns the ghost node’s index given its
1171    connectivity entry. In Fortran, use the expression :math:`-i`.
1172
1173 -  FEM_To_ghost_index(i) returns the connectivity entry for a given
1174    ghost node index. In Fortran, again use the expression :math:`-i`.
1175
1176 For example, a quadrilateral ghost element that is adjacent to,
1177 respectively, two real nodes 23 and 17, the tenth local ghost node, and
1178 one not-present node might have a connectivity entry of 23,17,-11,-1 (in
1179 C) or 23,17,-10,0 (in Fortran).
1180
1181 Applications may wish to use some other numbering, such as by storing
1182 all the ghost nodes after all the real nodes. The code to extract and
1183 renumber the connectivity of some 3-node triangles stored in FEM_ELEM+2
1184 would be:
1185
1186 ::
1187
1188    /* C version */
1189    int nReal=FEM_Mesh_get_length(mesh,FEM_ELEM+2);
1190    int nGhost=FEM_Mesh_get_length(mesh,FEM_GHOST+FEM_ELEM+2);
1191    typedef int intTriplet[3];
1192    intTriplet *conn=new intTriplet[nReal+nGhost];
1193    /* Extract real triangles into conn[0..nReal-1] */
1194    FEM_Mesh_data(mesh,FEM_ELEM+2,FEM_CONN, &conn[0][0], 0,nReal, 3,FEM_INDEX_0);
1195    /* Extract ghost triangles into conn[nReal..nReal+nGhost-1] */
1196    FEM_Mesh_data(mesh,FEM_GHOST+FEM_ELEM+2,FEM_CONN, &conn[nReal][0], 0,nGhost, 3,FEM_INDEX_0);
1197
1198    /* Renumber the ghost triangle connectivity */
1199    for (int t=nReal;t<nReal+nGhost;t++)
1200      for (int i=0;i<3;i++) {
1201        int in=conn[t][i]; /* uses FEM ghost node numbering */
1202        int out; /* uses application's ghost numbering */
1203        if (in==-1) {
1204          out=some_value_for_missing_nodes;
1205        } else if (FEM_Is_ghost_index(in)) {
1206          out=first_application_ghost+FEM_From_ghost_index(in);
1207        } else /*regular real node*/ {
1208          out=in;
1209        }
1210        conn[t][i]=out;
1211      }
1212
1213
1214 .. code-block:: fortran
1215
1216    ! F90 version
1217    INTEGER, ALLOCATABLE :: conn(3,:)
1218    INTEGER :: nReal,nGhost,t,i,in,out
1219    nReal=FEM_Mesh_get_length(mesh,FEM_ELEM+2)
1220    nGhost=FEM_Mesh_get_length(mesh,FEM_GHOST+FEM_ELEM+2)
1221    ALLOCATE(conn(3,nReal+nGhost))
1222    ! Extract real triangles into conn[1..nReal]
1223    CALL FEM_Mesh_data(mesh,FEM_ELEM+2,FEM_CONN, conn, 1,nReal, 3,FEM_INDEX_1)
1224    ! Extract ghost triangles into conn[nReal+1..nReal+nGhost]
1225    CALL FEM_Mesh_data(mesh,FEM_GHOST+FEM_ELEM+2,FEM_CONN, conn(1,nReal+1), 1,nGhost, 3,FEM_INDEX_1)
1226
1227    ! Renumber the ghost triangle connectivity
1228    DO t=nReal+1,nReal+nGhost
1229      DO i=1,3
1230        in=conn(i,t)
1231        IF (in .EQ. 0) out=some_value_for_missing_nodes
1232        IF (in .LT. 0) out=first_application_ghost-1+(-in)
1233        IF (in .GT. 0) out=in
1234        conn(i,t)=out
1235      END DO
1236    END DO
1237
1238
1239
1240 Setting up the ghost layer
1241 --------------------------
1242
1243 The framework’s ghost handling is element-centric. You specify which
1244 kinds of elements should be ghosts and how they connect by listing their
1245 faces before partitioning.
1246
1247 ::
1248
1249   void FEM_Add_ghost_layer(int nodesPerFace,int doAddNodes);
1250
1251 .. code-block:: fortran
1252
1253   SUBROUTINE FEM_Add_ghost_layer(nodesPerFace,doAddNodes)
1254   INTEGER, INTENT(IN) :: nodesPerFace,doAddNodes
1255
1256 This routine creates a new layer of ghosts around each FEM chunk.
1257 nodesPerFace is the number of shared nodes that together form a
1258 “face”. doAddNodes specifies that you want ghost nodes around your
1259 ghost elements. If doAddNodes is 0, ghost elements will have
1260 invalid -1 (in C) or 0 (in Fortran) connectivity entries where
1261 there is no corresponding local node.
1262
1263 A face is an unordered “tuple” of nodes, and is an abstract way to
1264 describe which ghosts your application needs—an element will be added
1265 to your chunk if it connects to at least one of your elements’ faces.
1266 For example, if you have a 3D, tetrahedral element that require
1267 ghosts on all 4 of its sides, this is equivalent to requiring ghosts
1268 of every element that shares all 3 nodes of one of your triangular
1269 faces, so for you a face is a 3-node triangle. If you have a 2D shape
1270 and want edge-adjacency, for you a face is a 2-node edge. If you want
1271 node-adjacent ghosts, a face is a single node.
1272
1273  Calling this routine several times creates several layers of ghost
1274  elements, and the different layers need not have the same parameters.
1275
1276 ::
1277
1278   void FEM_Add_ghost_elem(int elType,int facesPerElem,const int
1279      *elem2face);
1280
1281 .. code-block:: fortran
1282
1283   SUBROUTINE FEM_Add_ghost_elem(elType,facesPerElem,elem2face)
1284   INTEGER, INTENT(IN) :: elType,facesPerElem
1285   INTEGER, INTENT(IN) :: elem2face(nodesPerFace,facesPerElem)
1286
1287 This call is used to specify which type of element is to be added
1288 to the current ghost layer. facesPerElem and elem2face specify a
1289 mapping between each element and the surrounding faces. The
1290 elem2face table lists, for each face, the nodes of this element
1291 which form the face, specified as element-local numbers—indices
1292 into this element’s connectivity entry. The elem2face table should
1293 have nodesPerFace*facesPerElem entries, and no entry should be
1294 greater than nodePerEl for that element type.
1295
1296 Because all faces must take up the same space in the array, elem2face
1297 can include special indices— -1 for C, 0 for Fortran—that indicate
1298 the corresponding face is actually shorter than usual. For example,
1299 if nodesPerFace for this layer is 4, for 4-node quadrilateral faces,
1300 you could set one entry in elem2face to -1 to specify this is a
1301 3-node triangular face. Faces of different lengths will never match,
1302 so this is just a simple way to add ghosts from two kinds of faces at
1303 once.
1304
1305 The above two routines are always used together. For example, if your
1306 elements are 3-node triangles and you only require one shared node for
1307 inclusion in a single ghost layer, you would use:
1308
1309 ::
1310
1311   FEM_Add_ghost_layer(1,1); /* 1 node per face: node adjacency */
1312   const static int tri2node[]={0,1,2};
1313   FEM_Add_ghost_elem(0,3,tri2node); /* triangles are surrounded by 3 nodes */
1314
1315 If you require two shared nodes (a shared edge), the code will look
1316 like:
1317
1318 ::
1319
1320   FEM_Add_ghost_layer(2,1); /* 2 nodes per face: edge adjacency */
1321   const static int tri2edge[]={0,1,  1,2,  2,0};
1322   FEM_Add_ghost_elem(0,3,tri2edge); /*triangles are surrounded by 3 edges */
1323
1324 Symmetries and Ghosts-Geometric Layer
1325 -------------------------------------
1326
1327 The FEM framework can create ghosts not only of things that are on other
1328 processors, but also for various problem symmetries, like mirror
1329 reflection, and various types of periodicities. The interface for these
1330 ghosts is simple—you ask for the symmetries to be created, then you will
1331 get extra ghosts along each symmetry boundary. The symmetry ghosts are
1332 updated properly during any communication, even if the symmetry ghosts
1333 are ghosts of real local elements from the same chunk.
1334
1335 .. figure:: fig/sym_ghost.png
1336    :name: fig:symghost
1337    :width: 3in
1338
1339    Illustrating symmetry ghost elements.
1340
1341 Figure :numref:`fig:symghost` shows a chunk of a mesh for a rectangular
1342 domain with horizontal linear translational periodicity—that is, the
1343 domain repeats horizontally. Symmetry ghosts lie along the left and
1344 right sides; ordinary cross-processor parallel ghosts lie along the top
1345 edge where this chunk joins up with the rest of the domain; and the
1346 external boundary along the bottom of the chunk has no ghosts.
1347
1348 ::
1349
1350   void FEM_Add_linear_periodicity( int nFaces,int nPer, const int
1351   *facesA,const int *facesB, int nNodes,const double *nodeLocs );
1352
1353 .. code-block:: fortran
1354
1355   SUBROUTINE FEM_Add_linear_periodicity(nFaces,nPer,facesA,facesB,
1356   nNodes,nodeLocs)
1357   INTEGER, INTENT(IN) :: nFaces, nPer, nNodes
1358   INTEGER, INTENT(IN) :: facesA(nPer,nFaces), facesB(nPer,nFaces)
1359   double precision, INTENT(IN) :: nodeLocs(3,nNodes)
1360
1361 Make facesA and facesB match up under linear translation. Each face of
1362 facesA must match up with exactly one face of facesB, but both the
1363 faces and the nodes within a face can be permuted in any order—the
1364 order is recovered by matching 3d locations in the nodeLocs array.
1365
1366 This call can be repeated, for example if the domain is periodic along
1367 several directions. This call can only be issued from init().
1368
1369 ::
1370
1371   void FEM_Sym_coordinates(int elTypeOrMinusOne,double *locs);
1372
1373 .. code-block:: fortran
1374
1375   SUBROUTINE FEM_Sym_coordinates(elTypeOrZero,locs)
1376   INTEGER, INTENT(IN) :: elTypeOrZero
1377   double precision, intent(inout) :: locs(3,<number of items>)
1378
1379 This call adjusts the 3d locations listed in locs so they respect the
1380 symmetries of their corresponding item. If elTypeOrZero is an element
1381 type, the locations are adjusted to match with the corresponding
1382 element; if elTypeOrZero is zero, the locations are adjusted to match
1383 up with the corresponding node.
1384
1385 This call is needed because symmetry ghost nodes and elements initially
1386 have their original locations, which must be adjusted to respect the
1387 symmetry boundaries. Thus this call is needed both for initial location
1388 data (e.g., from FEM_Get_node_data) as well as any communicated location
1389 data (e.g., from FEM_Update_ghost_field).
1390
1391 This call can only be issued from driver().
1392
1393 Advanced Symmetries and Ghosts-Lower Layer
1394 ------------------------------------------
1395
1396 The geometric symmetry layer in the preceding section is actually a thin
1397 wrapper around this lower, more difficult to use layer.
1398
1399 ::
1400
1401   void FEM_Set_sym_nodes(const int *canon,const int *sym);
1402
1403 .. code-block:: fortran
1404
1405   SUBROUTINE FEM_Set_sym_nodes(canon,sym)
1406   INTEGER, INTENT(IN) :: canon(nNodes)
1407   INTEGER, INTENT(IN) :: sym(nNodes)
1408
1409 This call describes all possible symmetries in an extremely terse
1410 format. It can only be called from init(). The “canonicalization
1411 array” canon maps nodes to their canonical representative—if
1412 canon(\ :math:`i`)=canon(\ :math:`j`), nodes :math:`i` and :math:`j`
1413 are images of each other under some symmetry. The sym array has bits
1414 set for each symmetry boundary passing through a node.
1415
1416 For example, a 2d domain with 6 elements A, B, C, D, E, and F and 12
1417 nodes numbered 1-12 that is mirror-symmetric on the horizontal
1418 boundaries but periodic in the vertical boundaries would look like:
1419
1420 .. code-block:: none
1421
1422       D^'|  D^ |  E^ |  F^ |  F^`
1423       -  1  -  2  -  3  -  4  -
1424       A' |  A  |  B  |  C  |  C`
1425       -  5  -  6  -  7  -  8  -
1426       D' |  D  |  E  |  F  |  F`
1427       -  9  - 10  -  11 -  12 -
1428       Av'|  Av |  Bv |  Cv |  Cv`
1429
1430      v indicates the value has been shifted down (bottom boundary),
1431      ^ indicates the value has been shifted up (top boundary),
1432      ' indicates the value has been copied from the left (right boundary),
1433      ` indicates the value has been copied from the right (left boundary).
1434
1435 If we mark the left border with 1, the top with 2, the right with 4, and
1436 the bottom with 8, this situation is indicated by topologically pasting
1437 the top row to the bottom row by setting their canon entries equal, and
1438 marking each node with its symmetries.
1439
1440 ==== ===== =================
1441 Node canon sym
1442 ==== ===== =================
1443 1    1     3 (left + top)
1444 2    2     2 (top)
1445 3    3     2 (top)
1446 4    4     6 (top + right)
1447 5    5     1 (left)
1448 6    6     0 (none)
1449 7    7     0 (none)
1450 8    8     4 (right)
1451 9    1     9 (left+bottom)
1452 10   2     8 (bottom)
1453 11   3     8 (bottom)
1454 12   4     12 (bottom+right)
1455 ==== ===== =================
1456
1457 ::
1458
1459   void FEM_Get_sym(int elTypeOrMinusOne,int *destSym);
1460
1461
1462 .. code-block:: fortran
1463
1464   SUBROUTINE FEM_Get_sym(elTypeOrZero,destSym);
1465   INTEGER, INTENT(IN) :: elTypeOrMinusOne
1466   INTEGER, INTENT(OUT) :: destSym(nItems)
1467
1468 This call extracts the list of symmetry conditions that apply to an
1469 item type. If elType is an element type, it returns the symmetry
1470 conditions that apply to that element type; if elType is -1 (zero for
1471 Fortran), it returns the symmetry conditions that apply to the nodes.
1472 Symmetry conditions are normally only nonzero for ghost nodes and
1473 elements.
1474
1475 Mirror symmetry conditions are not yet supported, nor are multiple
1476 layers of symmetry ghosts, but both should be easy to add without
1477 changing this interface.
1478
1479 Older Mesh Routines
1480 ===================
1481
1482 These routines have a simpler, but less flexible interface than the
1483 general routines described in Section :numref:`sec:entities`. Because
1484 they are easy to implement in terms of the new routines, they will
1485 remain part of the framework indefinitely. These routines always use the
1486 default mesh, as returned by FEM_Mesh_default_read and
1487 FEM_Mesh_default_write.
1488
1489 ::
1490
1491   void FEM_Set_elem(int elType,int nEl,int doublePerEl,int nodePerEl);
1492   void FEM_Get_elem(int elType,int *nEl,int *doublePerEl,int
1493   *nodePerEl);
1494
1495 .. code-block:: fortran
1496
1497   SUBROUTINE FEM_Set_elem(elType,nEl,doublePerEl,nodePerEl)
1498   INTEGER, INTENT(IN) :: elType,nEl,doublePerEl,nodePerEl
1499   SUBROUTINE FEM_Get_elem(elType,nEl,doublePerEl,nodePerEl)
1500   INTEGER, INTENT(IN) :: elType
1501   INTEGER, INTENT(OUT) :: nEl,doublePerEl,nodePerEl
1502
1503 Describe/retrieve the number and type of elements. ElType is a
1504 user-defined small, unique element type tag. nEl is the number of
1505 elements being registered. doublesPerEl and nodePerEl are the number
1506 of doubles of user data, and nodes (respectively) associated with each
1507 element.
1508
1509 doublePerEl or nodePerEl may be zero, indicating that no user data or
1510 connectivity data (respectively) is associated with the element.
1511
1512 You can make this and any other mesh setup calls in any order—there is
1513 no need to make them in linearly increasing order. However, for a given
1514 type of element FEM_Set_elem must be called before setting that
1515 element’s connectivity or data.
1516
1517 ::
1518
1519   void FEM_Set_elem_conn(int elType,const int *conn);
1520   void FEM_Get_elem_conn(int elType,int *conn);
1521
1522 .. code-block:: fortran
1523
1524   SUBROUTINE FEM_Set_elem_conn_r(elType,conn)
1525   INTEGER, INTENT(IN) :: elType
1526   INTEGER, INTENT(IN), dimension(nodePerEl,nEl) :: conn
1527   SUBROUTINE FEM_Get_elem_conn_r(elType,conn)
1528   INTEGER, INTENT(IN) :: elType
1529   INTEGER, INTENT(OUT), dimension(nodePerEl,nEl) :: conn
1530   SUBROUTINE FEM_Set_elem_conn_c(elType,conn)
1531   INTEGER, INTENT(IN) :: elType
1532   INTEGER, INTENT(IN), dimension(nEl,nodePerEl) :: conn
1533   SUBROUTINE FEM_Get_elem_conn_c(elType,conn)
1534   INTEGER, INTENT(IN) :: elType
1535   INTEGER, INTENT(OUT), dimension(nEl,nodePerEl) :: conn
1536
1537 Describe/retrieve the element connectivity array for this element
1538 type. The connectivity array is indexed by the element number, and
1539 gives the indices of the nodes surrounding the element. It is hence
1540 nodePerEl*nEl integers long.
1541
1542 The C version array indices are zero-based, and must be stored in
1543 row-major order (a given element’s surrounding nodes are stored
1544 contiguously in the conn array). The Fortran version indices are
1545 one-based, and are available in row-major (named \_r) and column-major
1546 (named \_c) versions. We recommend row-major storage because it results
1547 in better cache utilization (because the nodes around an element are
1548 stored contiguously).
1549
1550 In this older interface, ghost nodes are indicated by invalid,
1551
1552 ::
1553
1554   void FEM_Set_node(int nNode,int doublePerNode);
1555   void FEM_Get_node(int
1556     *nNode,int *doublePerNode);
1557
1558 .. code-block:: fortran
1559
1560   SUBROUTINE FEM_Set_node(nNode,doublePerNode)
1561   INTEGER, INTENT(IN) :: nNode,doublePerNode
1562   SUBROUTINE FEM_Get_node(nNode,doublePerNode)
1563   INTEGER, INTENT(OUT) :: nNode,doublePerNode
1564
1565 Describe/retrieve the number of nodes and doubles of user data
1566 associated with each node. There is only one type of node, so no
1567 nodeType identifier is needed.
1568
1569 doublePerNode may be zero, indicating that no user data is associated
1570 with each node.
1571
1572 Old Mesh Data
1573 -------------
1574
1575 ::
1576
1577   void FEM_Set_node_data(const double *data);
1578   void FEM_Get_node_data(double *data);
1579   void FEM_Set_elem_data(int elType,const double *data);
1580   void FEM_Get_elem_data(int elType,double *data);
1581
1582 .. code-block:: fortran
1583
1584   SUBROUTINE FEM_Set_node_data_r(data)
1585   REAL*8, INTENT(IN), dimension(doublePerNode,nNode) :: data
1586   SUBROUTINE FEM_Get_node_data_r(data)
1587   REAL*8, INTENT(OUT), dimension(doublePerNode,nNode) :: data
1588   SUBROUTINE FEM_Set_node_data_c(data)
1589   REAL*8, INTENT(IN), dimension(nNode,doublePerNode) :: data
1590   SUBROUTINE FEM_Get_node_data_c(data)
1591   REAL*8, INTENT(OUT), dimension(nNode,doublePerNode) :: data
1592   SUBROUTINE FEM_Set_elem_data_r(elType,data)
1593   INTEGER, INTENT(IN) :: elType
1594   REAL*8, INTENT(IN), dimension(doublePerElem,nElem) :: data
1595   SUBROUTINE FEM_Get_elem_data_r(elType,data)
1596   INTEGER, INTENT(IN) :: elType
1597   REAL*8, INTENT(OUT), dimension(doublePerElem,nElem) :: data
1598   SUBROUTINE FEM_Set_elem_data_c(elType,data)
1599   INTEGER, INTENT(IN) :: elType
1600   REAL*8, INTENT(IN), dimension(nElem,doublePerElem) :: data
1601   SUBROUTINE FEM_Get_elem_data_c(elType,data)
1602   INTEGER, INTENT(IN) :: elType
1603   REAL*8, INTENT(OUT), dimension(nElem,doublePerElem) :: data
1604
1605 Describe/retrieve the optional, uninterpreted user data associated
1606 with each node and element. This user data is partitioned and
1607 reassembled along with the connectivity matrix, and may include
1608 initial conditions, node locations, material types, or any other data
1609 needed or produced by the program. The Fortran arrays can be row- or
1610 column- major (see FEM_Set_elem_conn for details). The row-major form
1611 is preferred.
1612
1613 Old Ghost Numbering
1614 -------------------
1615
1616 In this older version of the framework, FEM_Get_node and FEM_Get_elem
1617 return the **total** number of nodes and elements, including ghosts. The
1618 routines below return the index of the first ghost node or element,
1619 where ghosts are numbered after all the real elements. This old ghost
1620 numbering scheme does not work well when adding new ghosts, which is why
1621 the new ghost numbering scheme describes in
1622 Section :numref:`sec:ghostnum` is used in the new API.
1623
1624 .. figure:: fig/conn_indexing_old.png
1625    :name: fig:connold
1626    :width: 4in
1627
1628    Old ghost element and node numbering. FEM_Get_ghost_returns
1629    :math:`g`, FEM_Get_returns :math:`n`.
1630
1631 ::
1632
1633   int FEM_Get_node_ghost(void);
1634   int FEM_Get_elem_ghost(int elemType);
1635
1636 The examples below iterate over the real and ghost elements using the
1637 old numbering:
1638
1639 ::
1640
1641   // C version:
1642   int firstGhost,max;
1643   FEM_Get_node(&max, &ignored);
1644   firstGhost=FEM_Get_node_ghost();
1645   for (i=0;i<firstGhost;i++)
1646          //... i is a real node...
1647   for (i=firstGhost;i<max;i++)
1648          //... i is a ghost node ...
1649
1650 .. code-block:: fortran
1651
1652   ! Fortran version:
1653   call FEM_Get_node(max,ignored);
1654   firstGhost=FEM_Get_node_ghost();
1655   do i=1,firstGhost-1
1656   !       ... i is a real node...
1657   end do
1658   do i=firstGhost,max
1659   !      ... i is a ghost node...
1660   end do
1661
1662 Old Backward Compatibility
1663 --------------------------
1664
1665 ::
1666
1667   void FEM_Set_mesh(int nElem, int nNodes, int nodePerEl,const int*
1668   conn);
1669
1670 This is a convenience routine equivalent to:
1671
1672 ::
1673
1674   FEM_Set_node(nNodes,0);
1675   FEM_Set_elem(0,nElem,0,nodePerEl);
1676   FEM_Set_elem_Conn(0,conn);
1677
1678 .. code-block:: fortran
1679
1680   SUBROUTINE FEM_Set_mesh(nElem,nNodes,nodePerEl,conn)
1681
1682 This is a convenience routine equivalent to:
1683
1684 .. code-block:: fortran
1685
1686   CALL FEM_Set_node(nNodes,0)
1687   CALL FEM_Set_elem(1,nElem,0,nodePerEl)
1688   CALL FEM_Set_elem_Conn_c(1,conn)
1689
1690 Old Sparse Data
1691 ---------------
1692
1693 Sparse data is typically used to represent boundary conditions. For
1694 example, in a structural dynamics program typically some nodes have an
1695 imposed force or position. The routines in this section are used to
1696 describe this kind of mesh-associated data—data that only applies to
1697 some “sparse” subset of the nodes or elements.
1698
1699 ::
1700
1701   void FEM_Set_sparse(int S_id,int nRec, const int *nodes,int
1702   nodesPerRec, const void *data,int dataPerRec,int dataType);
1703
1704 .. code-block:: fortran
1705
1706   SUBROUTINE FEM_Set_sparse(S_id,nRec,nodes,nodesPerRec,data,dataPerRec,dataType)
1707   INTEGER, INTENT(IN) :: S_id,nRec,nodesPerRec,dataPerRec,dataType
1708   INTEGER, INTENT(IN) :: nodes(nodesPerRec,nRec)
1709   varies, INTENT(IN) :: data(dataPerRec,nRec)
1710
1711 Register nRec sparse data records with the framework under the number
1712 S_id. The first call to FEM_Set_sparse must give a S_id of zero in C
1713 (1 in fortran); and subsequent calls to FEM_Set_sparse must give
1714 increasing consecutive S_ids.
1715
1716 One sparse data record consists of some number of nodes, listed in the
1717 nodes array, and some amount of user data, listed in the data array.
1718 Sparse data records are copied into the chunks that contains all that
1719 record’s listed nodes. Sparse data records are normally used to describe
1720 mesh boundary conditions- for node-associated boundary conditions,
1721 nodesPerRec is 1; for triangle-associated boundary conditions,
1722 nodesPerRec is 3.
1723
1724 In general, nodePerRec gives the number of nodes associated with each
1725 sparse data record, and nodes gives the actual node numbers. dataPerRec
1726 gives the number of data items associated with each sparse data record,
1727 and dataType, one of FEM_BYTE, FEM_INT, FEM_REAL, or FEM_DOUBLE, gives
1728 the type of each data item. As usual, you may change or delete the nodes
1729 and data arrays after this call returns.
1730
1731 For example, if the first set of sparse data is 17 sparse data records,
1732 each containing 2 nodes stored in bNodes and 3 integers stored in bDesc,
1733 we would make the call:
1734
1735 ::
1736
1737   /*C version*/
1738   FEM_Set_sparse(0,17, bNodes,2, bDesc,3,FEM_INT);
1739
1740 .. code-block:: fortran
1741
1742   ! Fortran version
1743   CALL FEM_Set_sparse(1,17, bNodes,2, bDesc,3,FEM_INT)
1744
1745 ::
1746
1747   void FEM_Set_sparse_elem(int S_id,const int *rec2elem);
1748
1749 .. code-block:: fortran
1750
1751   SUBROUTINE FEM_Set_sparse_elem(S_id,rec2elem)
1752   INTEGER, INTENT(IN) :: S_id
1753   INTEGER, INTENT(IN) :: rec2elem(2,nRec)
1754
1755 Attach the previously-set sparse records S_id to the given elements.
1756 rec2elem consists of pairs of integers—one for each sparse data
1757 record. The first integer in the pair is the element type to attach
1758 the sparse record to, and the second integer gives the element number
1759 within that type. For example, to attach the 3 sparse records at S_id
1760 to the elements numbered 10, 11, and 12 of the element type elType,
1761 use:
1762
1763 ::
1764
1765   /*C version*/
1766   int rec2elem[]={elType,10, elType,11, elType,12};
1767   FEM_Set_sparse_elem(S_id,rec2elem);
1768
1769 .. code-block:: fortran
1770
1771   ! Fortran version
1772   integer :: rec2elem(2,3);
1773   rec2elem(1,:)=elType
1774   rec2elem(2,1)=10; rec2elem(2,2)=11; rec2elem(2,3)=12;
1775   CALL FEM_Set_sparse_elem(S_id,rec2elem)
1776
1777 ::
1778
1779   int FEM_Get_sparse_length(int S_id);
1780   void FEM_Get_sparse(int S_id,int *nodes,void *data);
1781
1782 .. code-block:: fortran
1783
1784   function FEM_Get_sparse_length(S_id);
1785   INTEGER, INTENT(IN) :: S_id
1786   INTEGER, INTENT(OUT) :: FEM_Get_sparse_Length
1787   SUBROUTINE FEM_Get_sparse(S_id,nodes,data);
1788   INTEGER, INTENT(IN) :: S_id
1789   INTEGER, INTENT(OUT) :: nodes(nodesPerRec,FEM_Get_sparse_Length(S_id))
1790   varies, INTENT(OUT) :: data(dataPerRec,FEM_Get_sparse_Length(S_id))
1791
1792 Retrieve the previously registered sparse data from the framework.
1793 FEM_Get_sparse_length returns the number of records of sparse data
1794 registered under the given S_id; zero indicates no records are
1795 available. FEM_Get_sparse returns you the actual nodes (translated to
1796 local node numbers) and unchanged user data for these sparse records.
1797
1798 In this old interface, there is no way to access sparse ghosts.
1799
1800 Mesh Modification
1801 =================
1802
1803 ::
1804
1805   void FEM_Update_mesh(FEM_Update_mesh_fn routine, int
1806   callMeshUpdated,int doWhat);
1807
1808 .. code-block:: fortran
1809
1810   SUBROUTINE FEM_Update_mesh(routine,callMeshUpdated,doWhat)
1811   external, INTENT(IN) :: routine
1812   INTEGER, INTENT(IN) :: callMeshUpdated,doWhat
1813
1814 Reassemble the mesh chunks from each partition into a single serial
1815 mesh, and call the given routine on the assembled mesh. In this
1816 routine, which runs on processor 0, the FEM_Get and FEM_Set routines
1817 can manipulate the serial mesh. The parameter callMeshUpdated, which
1818 must be non-zero, is passed down to routine as
1819 routine(callMeshUpdated).
1820
1821 FEM_Get calls from driver() will only return the new mesh after a
1822 FEM_Update_mesh call where doWhat is FEM_MESH_UPDATE; otherwise FEM_Get
1823 from driver() will still return the old mesh. FEM_Update_mesh can only
1824 be called from driver; and must be called by the driver routine for
1825 every chunk.
1826
1827 ================= ======= ============ ======================================
1828 doWhat            Numeric Repartition? FEM_Update_mesh
1829 ================= ======= ============ ======================================
1830 FEM_MESH_OUTPUT   0       No           driver() continues alongside routine
1831 FEM_MESH_FINALIZE 2       No           driver() blocks until routine finishes
1832 FEM_MESH_UPDATE   1       Yes          driver() blocks for the new partition
1833 ================= ======= ============ ======================================
1834
1835 For example, FEM_Update_mesh(my_output_routine, k, FEM_MESH_OUTPUT)
1836 reassembles the mesh and calls a routine named my_output_routine(k)
1837 while the driver routines continue with the computation. This might be
1838 useful, for example, for writing out intermediate solutions as a single
1839 file; writing outputs from driver() is more efficient but often results
1840 in a separate file for each mesh chunk.
1841
1842 To block the driver routines during a call to a routine named
1843 my_finalize_routine(k), such as at the end of the computation when the
1844 drivers have no other work to do, use
1845 FEM_Update_mesh(my_finalize_routine, k, FEM_MESH_FINALIZE).
1846
1847 To reassemble, modify, and repartition the mesh, use
1848 FEM_Update_mesh(my_update_routine, k, FEM_MESH_UPDATE). It may be easier
1849 to perform major mesh modifications from my_update_routine(k) than the
1850 drivers, since the entire serial mesh is available to
1851 my_update_routine(k).
1852
1853 FEM_Update_mesh reassembles the serial mesh with an attempt to preserve
1854 the element and node global numbering. If the new mesh has the same
1855 number and type of elements and nodes, the global numbers (and hence
1856 serial mesh) will be unchanged. If new elements or nodes are added at
1857 each chunk, they will be assigned new unique global numbers. If elements
1858 or nodes are removed, their global numbers are not re-used- you can
1859 detect the resulting holes in the serial mesh since the user data
1860 associated with the deleted elements will be all zero. Generally,
1861 however, it is less error-prone to perform mesh modifications only in
1862 driver() or only in an update routine, rather than some in both.
1863
1864 IDXL Communication
1865 ==================
1866
1867 The FEM framework’s communication layer is called IDXL. This small
1868 library handles sending and receiving data to and from a sparse subset
1869 of 1D indices into a user array. The sparse index subset is called an
1870 "Index List", hence the name of the library.
1871
1872 .. _sec:IDXL:
1873
1874 Index Lists
1875 -----------
1876
1877 An Index List is the fundamental data structure of the IDXL library—for
1878 example, the list of shared nodes is an Index List. IDXL includes
1879 routines for building, combining, and sending and receiving Index Lists.
1880
1881 An Index List, as you might expect, is a list of indices that need to be
1882 sent and received. An Index List includes both the indices that need to
1883 be sent, as well as the indices to be received, from each chunk.
1884
1885 Consider two chunks :math:`a` and :math:`b` where :math:`b` needs some
1886 information :math:`a` has, such as if :math:`b` has ghosts of real
1887 elements on :math:`a`. :math:`a`\ ’s Index List thus has a send portion
1888 with the :math:`a`-local indices for the elements :math:`a` sends; and
1889 :math:`b`\ ’s Index List contains a receive portion with the
1890 :math:`b`-local indices for the elements :math:`b` receives. Thus across
1891 processors, the corresponding send and receive portions of :math:`a` and
1892 :math:`b`\ ’s Index Lists match, as shown in
1893 Figure :numref:`fig:indexlists`.
1894
1895 .. figure:: fig/indexlists.png
1896    :name: fig:indexlists
1897    :width: 5in
1898
1899    Illustrating how Index Lists match up :math:`a`\ ’s source elements
1900    with :math:`b`\ ’s ghost elements.
1901
1902 Index List Calls
1903 ~~~~~~~~~~~~~~~~
1904
1905 You refer to an Index List via an opaque handle—in C, the integer
1906 typedef IDXL_t; in Fortran, a bare INTEGER.
1907
1908 ::
1909
1910   IDXL_t FEM_Comm_shared(int mesh,int entity);
1911
1912 .. code-block:: fortran
1913
1914   INTEGER function FEM_Comm_shared(mesh,entity)
1915   INTEGER, INTENT(IN) :: mesh,entity
1916
1917 Return a read-only copy of the Index List of shared nodes. The send
1918 and receive portions of this list are identical, because each shared
1919 node is both sent and received. Shared nodes are most often used with
1920 the send/sum communication pattern.
1921
1922 Must be called from driver. mesh must be a reading mesh. entity must be
1923 FEM_NODE. You may not call IDXL_Destroy on the returned list.
1924
1925 ::
1926
1927   IDXL_t FEM_Comm_ghost(int mesh,int entity);
1928
1929 .. code-block:: fortran
1930
1931   INTEGER function FEM_Comm_ghost(mesh,entity)
1932   INTEGER, INTENT(IN) :: mesh,entity
1933
1934 Return a read-only copy of the Index List of ghost entities. The send
1935 portion of this list contains real, interior entities, which are sent
1936 away; the receive portion of the list contains the ghost entites,
1937 which are received. Ghosts are most often used with the send/recv
1938 communication pattern.
1939
1940 Elements to be sent out are listed starting at zero (one in Fortran);
1941 but ghost elements to be received are also listed starting at zero (one
1942 in Fortran). If real and ghost elements are kept in separate arrays,
1943 this is usable as-is; but if ghosts and real elements are kept together,
1944 you will need to shift the ghost indices using IDXL_Combine or
1945 IDXL_Shift.
1946
1947 This routine must be called from driver. mesh must be a reading mesh.
1948 entity must not include FEM_GHOST-ghosts are already included. You may
1949 not call IDXL_Destroy on the returned list.
1950
1951 ::
1952
1953   IDXL_t IDXL_Create(void);
1954
1955 .. code-block:: fortran
1956
1957   INTEGER function IDXL_Create()
1958
1959 Create a new, empty Index List. This list can then be filled up using
1960 IDXL_Copy or IDXL_Combine.
1961
1962 Must be called from driver. You must eventually call IDXL_Destroy on the
1963 returned list.
1964
1965 ::
1966
1967   void IDXL_Combine(IDXL_t dest,IDXL_t src,int startSend,int startRecv);
1968
1969 .. code-block:: fortran
1970
1971   SUBROUTINE IDXL_Combine(dest,src,startSend,startRecv)
1972   INTEGER, INTENT(IN) :: dest,src,startSend,startRecv
1973
1974 Add the shifted contents of the src Index List to dest. The send
1975 portion of src is shifted so the first index sent will be startSend;
1976 for a ghost index list this is the index of the first sent real
1977 entity. The receive portion of src is similarly shifted so the first
1978 index received will be startRecv; for a ghost index list this is the
1979 index of the first received ghost entity.
1980
1981 This routine does not check for duplicates—if an index originally
1982 appears in dest and the also in the shifted src, it will be listed
1983 twice.
1984
1985 Advanced Index List Calls
1986 ~~~~~~~~~~~~~~~~~~~~~~~~~
1987
1988 ::
1989
1990   void IDXL_Destroy(IDXL_t l);
1991
1992 .. code-block:: fortran
1993
1994   SUBROUTINE IDXL_Destroy(l)
1995   INTEGER, INTENT(IN) :: l
1996
1997 Destroy this Index List, and free the list storage allocated by the
1998 framework. Only call this routine with lists you created using
1999 IDXL_Create; not lists obtained directly from the FEM framework.
2000
2001 ::
2002
2003   void IDXL_Print(IDXL_t l);
2004
2005 .. code-block:: fortran
2006
2007   SUBROUTINE IDXL_Print(l)
2008   INTEGER, INTENT(IN) :: l
2009
2010 Print out the contents of this Index List. This routine shows both the
2011 send and receive indices on the list, for each chunk we communicate
2012 with.
2013
2014 ::
2015
2016   void IDXL_Copy(IDXL_t dest,IDXL_t src);
2017
2018 .. code-block:: fortran
2019
2020   SUBROUTINE IDXL_Print(dest,src)
2021   INTEGER, INTENT(IN) :: dest,src
2022
2023 Copy the contents of the source Index List into the destination Index
2024 List, which should be empty.
2025
2026 ::
2027
2028   void IDXL_Shift(IDXL_t l,int startSend,int startRecv);
2029
2030 .. code-block:: fortran
2031
2032   SUBROUTINE IDXL_Shift(l,startSend,startRecv)
2033   INTEGER, INTENT(IN) :: l,startSend,startRecv
2034
2035 Like IDXL_Combine, but only shifts the indices within a single list.
2036
2037 ::
2038
2039   void IDXL_Add_entity(int newIdx,int nBetween,int *between);
2040
2041 .. code-block:: fortran
2042
2043   SUBROUTINE IDXL_Add_node(newIdx,nBetween,between)
2044   INTEGER, INTENT(IN) :: newIdx,nBetween
2045   INTEGER, INTENT(IN) :: between(nBetween)
2046
2047 This call adds a new entity, with local index newIdx, to this Index
2048 List. The new entity is sent or received by each chunk that sends or
2049 receives all the entities listed in the between array. For example,
2050 when adding a new node along an edge, nBetween is 2 and between lists
2051 the endpoints of the edge; this way if the edge is shared with some
2052 chunk, the new node will be shared with that chunk.
2053
2054 This routine only affects the current chunk- no other chunks are
2055 affected. To ensure the communication lists match, IDXL_Add_entity must
2056 be called on all the chunks that send or receive the entity, to create
2057 the local copies of the entity.
2058
2059 IDXL_Add_entity adds the new entity to the end of the communication
2060 list, and so must be called in the same order on all the chunks that
2061 share the new entity. For example, if two new nodes :math:`x` and
2062 :math:`y` are added between chunks :math:`a` and :math:`b`, if chunk
2063 :math:`a` calls IDXL_Add_entity with its local number for :math:`x`
2064 before it calls IDXL_Add_entity with its local number for :math:`y`,
2065 chunk :math:`b` must also add its copy of node :math:`x` before adding
2066 :math:`y`.
2067
2068 .. _sec:IDXLLayout:
2069
2070 Data Layout
2071 -----------
2072
2073 IDXL is designed to send and receive data directly out of your arrays,
2074 with no intermediate copying. This means IDXL needs a completely general
2075 method for specifying how you store your data in your arrays. Since you
2076 probably don’t change your storage layout at runtime, you can create a
2077 “data layout” once at the beginning of your program, then use it
2078 repeatedly for communication.
2079
2080 IDXL Layouts are normally used to describe arrays of data associated
2081 with nodes or elements. The layout abstraction allows you to use IDXL
2082 routines to communicate any sort of data, stored in a variety of
2083 formats.
2084
2085 Like Index Lists, Layouts are referred to via an opaque handle—in a C
2086 program via the integer typedef IDXL_Layout_t, and in Fortran via a bare
2087 integer.
2088
2089 Layout Routines
2090 ~~~~~~~~~~~~~~~
2091
2092 In most programs, the data to be communicated is a dense array of data
2093 of one type. In this case, there is only one layout routine you need to
2094 know:
2095
2096 ::
2097
2098   IDXL_Layout_t IDXL_Layout_create(int type,int width);
2099
2100 .. code-block:: fortran
2101
2102   INTEGER function IDXL_Layout_create(type,width)
2103   INTEGER, INTENT(IN) :: type,width
2104
2105 The simplest data layout to describe—a dense array of this IDXL
2106 datatype, indexed by entity number, with width pieces of data per
2107 entity. Note that the number of entities is not stored with the
2108 layout-the number of entities to be communicated depends on the
2109 communication routine.
2110
2111 The IDXL datatypes are:
2112
2113 ============= ============= =================
2114 IDXL Datatype C Datatypes   Fortran Datatypes
2115 ============= ============= =================
2116 IDXL_BYTE     unsigned char INTEGER*1
2117 \             char          LOGICAL*1
2118 IDXL_INT      int           INTEGER
2119 IDXL_REAL     float         SINGLE PRECISION
2120 \                           REAL*4
2121 IDXL_DOUBLE   double        DOUBLE PRECISION
2122 \                           REAL*8
2123 ============= ============= =================
2124
2125 For example, if you keep a dense array with 3 doubles of force per node,
2126 you’d call this routine as:
2127
2128 ::
2129
2130   // C++ version:
2131   double *force=new double[3*n];
2132   IDXL_Layout_t force_layout=IDXL_Layout_create(IDXL_DOUBLE,3);
2133
2134 .. code-block:: fortran
2135
2136   ! F90 Version
2137   double precision, allocatable :: force(:,:)
2138   integer :: force_layout
2139   ALLOCATE(force(3,n)) ! (could equivalently use force(3*n) )
2140   force_layout=IDXL_Layout_create(IDXL_DOUBLE,3)
2141
2142 This routine was once called FEM_Create_simple_field.
2143
2144 .. _sec:IDXLLayoutoffset:
2145
2146 Advanced Layout Routines
2147 ~~~~~~~~~~~~~~~~~~~~~~~~
2148
2149 These advanced routines are only needed if you want to exchange data
2150 stored in an array of user-defined types. Most programs only need
2151 IDXL_Layout_create.
2152
2153 ::
2154
2155   IDXL_Layout_t IDXL_Layout_offset(int type, int width, int offsetBytes,
2156   int distanceBytes,int skewBytes);
2157
2158 .. code-block:: fortran
2159
2160   INTEGER function IDXL_Layout_offset(type,width,offsetBytes,distanceBytes,skewBytes)
2161   INTEGER, INTENT(IN) :: type,width,offsetBytes,distanceBytes,skewBytes
2162
2163 The most general data layout-an array indexed by entity, containing
2164 width pieces of data per entity. This routine expands on
2165 IDXL_Layout_create by adding support for user-defined types or other
2166 unusual data layouts. You describe your layout by giving various
2167 in-memory byte offsets that describe the data is stored. Again, the
2168 number of entities is not stored with the layout-the number of
2169 entities to be communicated depends on the communication routine.
2170
2171 -  offsetBytes The number of bytes from the start of the array to the
2172    start of the data.
2173
2174 -  distanceBytes The number of bytes taken by one entity.
2175
2176 -  skewBytes The number of bytes between each piece of data. Since this
2177    can almost always be determined from the size of the base data type,
2178    this parameter can be left as zero.
2179
2180 .. figure:: fig/layout.png
2181    :name: fig:layout
2182    :width: 5in
2183
2184    Describing a complex data layout.
2185
2186 For example, if your node data is all stored in a struct (in fortran, a
2187 named TYPE), offsetBytes gives the distance between the start of the
2188 struct and the force; and distanceBytes gives the size in bytes of the
2189 struct.
2190
2191 In C, the offsetof and sizeof keywords are useful for finding these
2192 values. In Fortran, we provide a special routine called foffsetof that
2193 returns the distance, in bytes, between its two arguments.
2194
2195 ::
2196
2197   // C++ version:
2198   typedef struct {
2199      double d[3], v[3], force[3], a[3];
2200      double m;
2201   } node;
2202   node *nodes=new node[n];
2203   IDXL_Layout_t force_layout=IDXL_Layout_offset(IDXL_DOUBLE,3,
2204            offsetof(node,force),sizeof(node),0);
2205
2206
2207 .. code-block:: fortran
2208
2209   ! F90 Version
2210     TYPE node
2211        DOUBLE PRECISION :: d(3), v(3), force(3), a(3)
2212        DOUBLE PRECISION :: m
2213     END TYPE
2214     integer :: force_layout
2215     ALLOCATE(nodes(n))
2216     force_layout=IDXL_Layout_create(IDXL_DOUBLE,3,
2217   &          foffsetof(nodes(1),nodes(1)%force),
2218   &          foffsetof(nodes(1),nodes(2)),0)
2219
2220 ::
2221
2222   void IDXL_Layout_destroy(IDXL_Layout_t layout);
2223
2224 .. code-block:: fortran
2225
2226   SUBROUTINE IDXL_Layout_destroy(layout)
2227   INTEGER, INTENT(IN) :: layout
2228
2229 Destroy this Layout. You only need call this routine if you repeatedly
2230 create layouts.
2231
2232 ::
2233
2234   int IDXL_Get_layout_type(IDXL_Layout_t layout);
2235
2236 .. code-block:: fortran
2237
2238   INTEGER function IDXL_Get_layout_type(layout)
2239
2240 Return the IDXL datatype for this layout.
2241
2242 ::
2243
2244   int IDXL_Get_layout_width(IDXL_Layout_t layout);
2245
2246 .. code-block:: fortran
2247
2248   INTEGER function IDXL_Get_layout_width(layout)
2249
2250 Return the layout width—the number of data items that are communicated
2251 per entity.
2252
2253 ::
2254
2255   int IDXL_Get_layout_distance(IDXL_Layout_t layout);
2256
2257 .. code-block:: fortran
2258
2259   INTEGER function IDXL_Get_layout_distance(layout)
2260
2261 Return the layout distance—the number of bytes between successive
2262 entity’s data items.
2263
2264 Layout Compatibility Routines
2265 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2266
2267 Before IDXL was made a separate library, FEM included these routines,
2268 which are still preserved for backward compatibility.
2269
2270 ::
2271
2272   IDXL_Layout_t FEM_Create_simple_field(int type,int width);
2273
2274 .. code-block:: fortran
2275
2276   INTEGER function FEM_Create_simple_field(type,width)
2277   INTEGER, INTENT(IN) :: type,width
2278
2279 This routine is completely interchangeable to IDXL_Layout_create.
2280
2281 ::
2282
2283   int FEM_Create_field(int type,int width,int offset,int distance);
2284
2285 .. code-block:: fortran
2286
2287   INTEGER function FEM_Create_field(type, width, offset, distance)
2288   INTEGER, INTENT(IN) :: type, width, offset, distance
2289
2290 This routine is like a call to IDXL_Layout_offset with the rarely used
2291 skewBytes set to zero.
2292
2293 .. _sec:IDXLComm:
2294
2295 IDXL Communication
2296 ------------------
2297
2298 This section brings together all the pieces of IDXL: Index Lists are
2299 used to determine what to send and what to receive and Layouts are used
2300 to determine where to get and put the communicated data.
2301
2302 Communication Routines
2303 ~~~~~~~~~~~~~~~~~~~~~~
2304
2305 ::
2306
2307   void IDXL_Comm_sendsum(IDXL_Comm_t comm,IDXL_t indices,IDXL_Layout_t
2308   layout,void *data);
2309
2310 .. code-block:: fortran
2311
2312   SUBROUTINE IDXL_Comm_sendsum(comm,indices,layout,data)
2313   INTEGER, INTENT(IN) :: comm,indices,layout
2314   varies, INTENT(INOUT) :: data
2315
2316 Sum these indices of shared entities across all chunks that share
2317 them. The user data array is interpreted according to the given
2318 layout.
2319
2320 If comm is zero, this routine is blocking and finishes the communication
2321 immediately. If comm is not zero, this routine is non-blocking and
2322 equivalent to a call to IDXL_Comm_send followed by a call to
2323 IDXL_Comm_sum.
2324
2325 This routine is typically used to sum up partial values on shared nodes.
2326 It is a more general version of the old FEM routine FEM_Update_field.
2327 For example, to sum up the shared-node values in a 3d force vector
2328 indexed by node, you would use:
2329
2330 ::
2331
2332   // C++ version:
2333   double *force=new double[3*nNodes];
2334   IDXL_Layout_t force_layout=IDXL_Layout_create(IDXL_DOUBLE,3);
2335   IDXL_t shared_indices=FEM_Comm_shared(mesh,FEM_NODE);
2336
2337   //... in the time loop ...
2338       IDXL_Comm_sendsum(0,shared_indices,force_layout,force);
2339
2340 .. code-block:: fortran
2341
2342   ! F90 Version
2343   double precision, allocatable :: force(:,:)
2344   integer :: force_layout, shared_indices
2345   ALLOCATE(force(3,nNodes)) ! (could equivalently use force(3*nNodes) )
2346   force_layout=IDXL_Layout_create(IDXL_DOUBLE,3)
2347   shared_indices=FEM_Comm_shared(mesh,FEM_NODE)
2348
2349   !... in the time loop ...
2350       CALL IDXL_Comm_sendsum(0,shared_indices,force_layout,force)
2351
2352 ::
2353
2354   void IDXL_Comm_sendrecv(IDXL_Comm_t comm,IDXL_t indices,IDXL_Layout_t
2355   layout,void *data);
2356
2357 .. code-block:: fortran
2358
2359   SUBROUTINE IDXL_Comm_sendrecv(comm,indices,layout,data)
2360   INTEGER, INTENT(IN) :: comm,indices,layout
2361   varies, INTENT(INOUT) :: data
2362
2363 Send these (typically real) send indices and copy in these (typically
2364 ghost) receive indices. The user data array is interpreted according
2365 to the given layout.
2366
2367 If comm is zero, this routine is blocking and finishes the communication
2368 immediately. If comm is not zero, this routine is non-blocking and
2369 equivalent to a call to IDXL_Comm_send followed by a call to
2370 IDXL_Comm_sum.
2371
2372 This routine is typically used to obtain the values of ghost entities.
2373 It is a more general version of the old FEM routine
2374 FEM_Update_ghost_field. For example, to obtain 7 solution values per
2375 ghost element, storing gElem ghosts in the array just after the nElem
2376 regular elements, we could:
2377
2378 ::
2379
2380   // C++ version:
2381   double *elem=new double[7*(nElem+gElem)];
2382   IDXL_Layout_t elem_layout=IDXL_Layout_create(IDXL_DOUBLE,7);
2383   IDXL_t ghost_original=FEM_Comm_ghost(mesh,FEM_ELEM+1);
2384   IDXL_t ghost_shifted=IDXL_Create(); // ghosts start at nElem
2385   IDXL_Combine(ghost_shifted,ghost_original,0,nElem);
2386
2387   //... in the time loop ...
2388       IDXL_Comm_sendrecv(0,ghost_shifted,elem_layout,elem);
2389
2390 .. code-block:: fortran
2391
2392   ! F90 Version
2393   double precision, allocatable :: elem(:,:)
2394   integer :: elem_layout, ghost_original,ghost_shifted
2395   ALLOCATE(elem(7,nElem+gElem))
2396   elem_layout=IDXL_Layout_create(IDXL_DOUBLE,7)
2397   ghost_original=FEM_Comm_ghost(mesh,FEM_ELEM+1)
2398   ghost_shifted=IDXL_Create() ! ghosts start at nElem+1
2399   CALL IDXL_Combine(ghost_shifted,ghost_original,1,nElem+1)
2400
2401   !... in the time loop ...
2402       CALL IDXL_Comm_sendrecv(0,ghost_shifted,elem_layout,elem)
2403
2404 Advanced Communication Routines
2405 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2406
2407 ::
2408
2409   IDXL_Comm_t IDXL_Comm_begin(int tag,int context);
2410
2411 .. code-block:: fortran
2412
2413   INTEGER function IDXL_Comm_begin(tag,context)
2414   INTEGER, INTENT(IN) :: tag,context
2415
2416 Start a non-blocking communication operation with this (user-defined)
2417 tag and communication context (0, or an AMPI communicator).
2418
2419 Every call to this routine must eventually be matched by a call to
2420 IDXL_Comm_wait. Warning: for now, tag and context are ignored, and there
2421 can be only one outstanding communication operation.
2422
2423 ::
2424
2425   void IDXL_Comm_send(IDXL_Comm_t comm,IDXL_t indices,IDXL_Layout_t
2426   layout,const void *data);
2427
2428 .. code-block:: fortran
2429
2430   SUBROUTINE IDXL_Comm_send(comm,indices,layout,data)
2431   INTEGER, INTENT(IN) :: comm,indices,layout
2432   varies, INTENT(IN) :: data
2433
2434 When comm is flushed, send these send indices, with this layout, from
2435 this data array.
2436
2437 This routine is always non-blocking; as the data array passed in will
2438 not be copied out until the call to IDXL_Comm_flush.
2439
2440 ::
2441
2442   void IDXL_Comm_recv(IDXL_Comm_t comm,IDXL_t indices,IDXL_Layout_t
2443   layout,void *data);
2444
2445 .. code-block:: fortran
2446
2447   SUBROUTINE IDXL_Comm_recv(comm,indices,layout,data)
2448   INTEGER, INTENT(IN) :: comm,indices,layout
2449   varies, INTENT(OUT) :: data
2450
2451 When comm is finished, copy in these receive indices, with this
2452 layout, into this data array.
2453
2454 This routine is always non-blocking; as the data array passed in will
2455 not be copied into until the call to IDXL_Comm_wait.
2456
2457 ::
2458
2459   void IDXL_Comm_sum(IDXL_Comm_t comm,IDXL_t indices,IDXL_Layout_t
2460   layout,void *data);
2461
2462 .. code-block:: fortran
2463
2464   SUBROUTINE IDXL_Comm_sum(comm,indices,layout,data)
2465   INTEGER, INTENT(IN) :: comm,indices,layout
2466   varies, INTENT(INOUT) :: data
2467
2468 When comm is finished, add in the values for these receive indices,
2469 with this layout, into this data array.
2470
2471 This routine is always non-blocking; as the data array passed in will
2472 not be added to until the call to IDXL_Comm_wait.
2473
2474 ::
2475
2476   void IDXL_Comm_flush(IDXL_Comm_t comm);
2477
2478 .. code-block:: fortran
2479
2480   SUBROUTINE IDXL_Comm_flush(comm)
2481   INTEGER, INTENT(IN) :: comm
2482
2483 Send all outgoing data listed on this comm. This routine exists
2484 because there may be many calls to IDXL_Comm_send, and sending one
2485 large message is more efficient than sending many small messages.
2486
2487 This routine is typically non-blocking, and may only be issued at most
2488 once per IDXL_Comm_begin.
2489
2490 ::
2491
2492   void IDXL_Comm_wait(IDXL_Comm_t comm);
2493
2494 .. code-block:: fortran
2495
2496   SUBROUTINE IDXL_Comm_wait(comm)
2497   INTEGER, INTENT(IN) :: comm
2498
2499 Finish this communication operation. This call must be issued exactly
2500 once per IDXL_Comm_begin. This call includes IDXL_Comm_flush if it has
2501 not yet been called.
2502
2503 This routine always blocks until all incoming data is received, and is
2504 the last call that can be made on this comm.
2505
2506 Old Communication Routines
2507 ==========================
2508
2509 (This section is for backward compatibility only. The IDXL routines are
2510 the new, more flexible way to perform communication.)
2511
2512 The FEM framework handles the updating of the values of shared nodes-
2513 that is, it combines shared nodes’ values across all processors. The
2514 basic mechanism to do this update is the “field”- numeric data items
2515 associated with each node. We make no assumptions about the meaning of
2516 the node data, allow various data types, and allow a mix of communicated
2517 and non-communicated data associated with each node. The framework uses
2518 IDXL layouts to find the data items associated with each node in memory.
2519
2520 Each field represents a (set of) data records stored in a contiguous
2521 array, often indexed by node number. You create a field once, with the
2522 IDXL layout routines or FEM_Create_field, then pass the resulting field
2523 ID to FEM_Update_field (which does the shared node communication),
2524 FEM_Reduce_field (which applies a reduction over node values), or one of
2525 the other routines described below.
2526
2527 ::
2528
2529   void FEM_Update_field(int Fid,void *nodes);
2530
2531 .. code-block:: fortran
2532
2533   SUBROUTINE FEM_Update_field(Fid,nodes)
2534   INTEGER, INTENT(IN) :: Fid
2535   varies, INTENT(INOUT) :: nodes
2536
2537 Combine a field of all shared nodes with the other chunks. Sums the
2538 value of the given field across all chunks that share each node. For
2539 the example above, once each chunk has computed the net force on each
2540 local node, this routine will sum the net force across all shared
2541 nodes.
2542
2543 FEM_Update_field can only be called from driver, and to be useful, must
2544 be called from every chunk’s driver routine.
2545
2546 After this routine returns, the given field of each shared node will be
2547 the same across all processors that share the node.
2548
2549 This routine is equivalent to an IDXL_Comm_Sendsum operation.
2550
2551 ::
2552
2553   void FEM_Read_field(int Fid,void *nodes,char *fName);
2554
2555 .. code-block:: fortran
2556
2557   SUBROUTINE FEM_Read_field(Fid,nodes,fName)
2558   INTEGER, INTENT(IN) :: Fid
2559   varies, INTENT(OUT) :: nodes
2560   CHARACTER*, INTENT(IN) :: fName
2561
2562 Read a field out of the given serial input file. The serial input file
2563 is line-oriented ASCII- each line begins with the global node number
2564 (which must match the line order in the file), followed by the data to
2565 be read into the node field. The remainder of each line is unread. If
2566 called from Fortran, the first line must be numbered 1; if called from
2567 C, the first line must be numbered zero. All fields are separated by
2568 white space (any number of tabs or spaces).
2569
2570 For example, if we have called Create_field to describe 3 doubles, the
2571 input file could begin with
2572
2573 .. code-block:: none
2574
2575              1    0.2    0.7    -0.3      First node
2576              2    0.4    1.12   -17.26    another node
2577              ...
2578
2579 FEM_Read_field must be called from driver at any time, independent of
2580 other chunks.
2581
2582 This routine has no IDXL equivalent.
2583
2584 ::
2585
2586   void FEM_Reduce_field(int Fid,const void *nodes,void *out,int op);
2587
2588 .. code-block:: fortran
2589
2590   SUBROUTINE FEM_Reduce_field(Fid,nodes,outVal,op)
2591   INTEGER, INTENT(IN) :: Fid,op
2592   varies, INTENT(IN) :: nodes
2593   varies, INTENT(OUT) :: outVal
2594
2595 Combine one record per node of this field, according to op, across all
2596 chunks. Shared nodes are not double-counted- only one copy will
2597 contribute to the reduction. After Reduce_field returns, all chunks
2598 will have identical values in outVal, which must be vec_len copies of
2599 base_type.
2600
2601 May only be called from driver, and to complete, must be called from
2602 every chunk’s driver routine.
2603
2604 op must be one of:
2605
2606 -  FEM_SUM- each element of outVal will be the sum of the corresponding
2607    fields of all nodes
2608
2609 -  FEM_MIN- each element of outVal will be the smallest value among the
2610    corresponding field of all nodes
2611
2612 -  FEM_MAX- each element of outVal will be the largest value among the
2613    corresponding field of all nodes
2614
2615 This routine has no IDXL equivalent.
2616
2617 ::
2618
2619   void FEM_Reduce(int Fid,const void *inVal,void *outVal,int op);
2620
2621 .. code-block:: fortran
2622
2623   SUBROUTINE FEM_Reduce(Fid,inVal,outVal,op)
2624   INTEGER, INTENT(IN) :: Fid,op
2625   varies, INTENT(IN) :: inVal
2626   varies, INTENT(OUT) :: outVal
2627
2628 Combine one record of this field from each chunk, according to op,
2629 across all chunks. Fid is only used for the base_type and vec_len-
2630 offset and dist are not used. After this call returns, all chunks will
2631 have identical values in outVal. Op has the same values and meaning as
2632 FEM_Reduce_field.
2633
2634 May only be called from driver, and to complete, must be called from
2635 every chunk’s driver routine.
2636
2637 ::
2638
2639   ! C example
2640   double inArr[3], outArr[3];
2641   int fid=IDXL_Layout_create(FEM_DOUBLE,3);
2642   FEM_Reduce(fid,inArr,outArr,FEM_SUM);
2643
2644 .. code-block:: fortran
2645
2646   ! f90 example
2647   DOUBLE PRECISION :: inArr(3), outArr(3)
2648   INTEGER fid
2649   fid=IDXL_Layout_create(FEM_DOUBLE,3)
2650   CALL FEM_Reduce(fid,inArr,outArr,FEM_SUM)
2651
2652 This routine has no IDXL equivalent.
2653
2654 Ghost Communication
2655 -------------------
2656
2657 It is possible to get values for a chunk’s ghost nodes and elements from
2658 the neighbors. To do this, use:
2659
2660 ::
2661
2662   void FEM_Update_ghost_field(int Fid, int elTypeOrMinusOne, void
2663   *data);
2664
2665 .. code-block:: fortran
2666
2667   SUBROUTINE FEM_Update_ghost_field(Fid,elTypeOrZero,data)
2668   INTEGER, INTENT(IN) :: Fid,elTypeOrZero
2669   varies, INTENT(INOUT) :: data
2670
2671 This has the same requirements and call sequence as FEM_Update_field,
2672 except it applies to ghosts. You specify which type of element to
2673 exchange using the elType parameter. Specify -1 (C version) or 0
2674 (fortran version) to exchange node values.
2675
2676 Ghost List Exchange
2677 -------------------
2678
2679 It is possible to exchange sparse lists of ghost elements between FEM
2680 chunks.
2681
2682 ::
2683
2684   void FEM_Exchange_ghost_lists(int elemType,int nIdx,const int
2685   *localIdx);
2686
2687 .. code-block:: fortran
2688
2689   SUBROUTINE FEM_Exchange_ghost_lists(elemType,nIdx,localIdx)
2690   INTEGER, INTENT(IN) :: elemType,nIdx
2691   INTEGER, INTENT(IN) :: localIdx[nIdx]
2692
2693 This routine sends the local element indices in localIdx to those
2694 neighboring chunks that connect to its ghost elements on the other
2695 side. That is, if the element localIdx[i] has a ghost on some chunk c,
2696 localIdx[i] will be sent to and show up in the ghost list of chunk c.
2697
2698 ::
2699
2700   int FEM_Get_ghost_list_length(void);
2701
2702 Returns the number of entries in my
2703 ghost list—the number of my ghosts that other chunks passed to their
2704 call to FEM_Exchange_ghost_lists.
2705
2706 ::
2707
2708   void FEM_Get_ghost_list(int *retLocalIdx);
2709
2710 .. code-block:: fortran
2711
2712   SUBROUTINE FEM_Get_ghost_list(retLocalIdx)
2713   INTEGER, INTENT(OUT) :: retLocalIdx[FEM_Get_ghost_list_length()]
2714
2715 These routines access the list of local elements sent by other chunks.
2716   The returned indices will all refer to ghost elements in my chunk.
2717
2718 .. _sec:ParFUM:
2719
2720 ParFUM
2721 ======
2722
2723 ParFUM is the name for the latest version of FEM. ParFUM includes
2724 additional features including parallel mesh modification and adaptivity
2725 (geometrical). ParFUM also contains functions which generate additional
2726 topological adjacency information. ParFUM cannot be built separate from
2727 Charm++ since it uses various messaging mechanisms that MPI does not
2728 readily support. It is important to note that ParFUM adaptivity at the
2729 moment has some limitations. It works only for meshes which are
2730 two-dimensional. The other limitation is that the mesh on which it works
2731 on must have one layer of node-deep ghosts. Most applications require no
2732 or one layer ghosts, so it is really not a limitation, but for
2733 applications that need multiple layers of ghost information, the
2734 adaptivity operations cannot be used.
2735
2736 Adaptivity Initialization
2737 -------------------------
2738
2739 If a FEM application wants to use parallel mesh adaptivity, the first
2740 task is to call the initialization routine from the *driver* function.
2741 This creates the node and element adjacency information that is
2742 essential for the adaptivity operations. It also initializes all the
2743 mesh adaptivity related internal objects in the framework.
2744
2745 ::
2746
2747   void FEM_ADAPT_Init(int meshID)
2748
2749 Initializes the mesh defined by meshID for the mesh adaptivity
2750 operations.
2751
2752 Preparing the Mesh for Adaptivity
2753 ---------------------------------
2754
2755 For every element entity in the mesh, there is a desired size entry for
2756 each element. This entry is called meshSizing. This meshSizing entry
2757 contains a metric that decides the element quality. The default metric
2758 is the average of the size of the three edges of an element. This
2759 section provides various mechanisms to set this field. Some of the
2760 adaptive operations actually use these metrics to maintain quality.
2761 Though there is another metric which is computer for each element and
2762 maintained on the fly and that is the ratio of the largest length to the
2763 smallest altitude and this value during mesh adaptivity is not allowed
2764 to go beyond a certain limit. Because the larger this value after a
2765 certain limit, the worse the element quality.
2766
2767 ::
2768
2769   void FEM_ADAPT_SetElementSizeField(int meshID, int elem, double size);
2770
2771 For the mesh specified by meshID, for the element elem, we set the
2772 desired size for each element to be size.
2773
2774 ::
2775
2776   void FEM_ADAPT_SetElementSizeField(int meshID, double *sizes);
2777
2778 For the mesh specified by meshID, for the element elem, we set the
2779 desired size for each element from the corresponding entry in the sizes
2780 array.
2781
2782 ::
2783
2784   void FEM_ADAPT_SetReferenceMesh(int meshID);
2785
2786 For each element int this mesh defined by meshID set its size to the
2787 average edge length of the corresponding element.
2788
2789 ::
2790
2791   void FEM_ADAPT_GradateMesh(int meshID, double smoothness);
2792
2793 Resize mesh elements to avoid jumps in element size. i.e. avoid
2794 discontinuities in the desired sizes for elements of a mesh by smoothing
2795 them out. Algorithm based on h-shock correction, described in Mesh
2796 Gradation Control, Borouchaki et al.
2797
2798 Modifying the Mesh
2799 ------------------
2800
2801 Once the elements in the mesh has been prepared by specifying their
2802 desired sizes, we are ready to use the actual adaptivity operations.
2803 Currently we provide Delaunay flip operations, edge bisect operations
2804 and edge-coarsen operations all of which are implemented in parallel,
2805 but the user has access to these wrapper functions which intelligently
2806 decide when and in which region of the mesh to use the adaptivity
2807 operations to generate a mesh with higher quality elements while
2808 achieving the desired size (which is usually average edge length per
2809 element), or it could even be the area of each element.
2810
2811 ::
2812
2813   void FEM_ADAPT_Refine(int meshID, int qm, int method, double
2814   factor,double *sizes);
2815
2816 Perform refinements on the mesh specified by meshId. Tries to
2817 maintain/improve element quality by refining the mesh as specified by a
2818 quality measure qm. If method = 0, refine areas with size larger than
2819 factor down to factor If method = 1, refine elements down to sizes
2820 specified in the sizes array. In this array each entry corresponds to
2821 the corresponding element. Negative entries in sizes array indicate no
2822 refinement.
2823
2824 ::
2825
2826   void FEM_ADAPT_Coarsen(int meshID, int qm, int method, double
2827   factor,double *sizes);
2828
2829 Perform refinements on the mesh specified by meshId. Tries to
2830 maintain/improve element quality by coarsening the mesh as specified by
2831 a quality measure qm. If method = 0, coarsen areas with size smaller
2832 than factor down to factor If method = 1, coarsen elements up to sizes
2833 specified in the sizes array. In this array each entry corresponds to
2834 the corresponding element. Negative entries in sizes array indicate no
2835 coarsening.
2836
2837 ::
2838
2839   void FEM_ADAPT_AdaptMesh(int meshID, int qm, int method, double
2840   factor,double *sizes);
2841
2842 It has the same set of arguments as required by the previous two
2843 operations, namely refine and coarsen. This function keeps using the
2844 above two functions till we have all elements in the mesh with as close
2845 to the desired quality. Apart from using the above two operations, it
2846 also performs a mesh repair operation where it gets rid of some bad
2847 quality elements by Delaunay flip or coarsening as the geometry in the
2848 area demands.
2849
2850 ::
2851
2852   int FEM_ADAPT_SimpleRefineMesh(int meshID, double targetA, double xmin,
2853   double ymin, double xmax, double ymax);
2854
2855 A region is defined by (xmax, xmin, ymax, ymin) and the target area to
2856 be achieved for all elements in this region in the mesh specified by
2857 meshID is given by targetA. This function only performs a series of
2858 refinements on the elements in this region. If the area is larger, then
2859 no coarsening is done.
2860
2861 ::
2862
2863   int FEM_ADAPT_SimpleCoarsenMesh(int meshID, double targetA, double xmin,
2864   double ymin, double xmax, double ymax);
2865
2866 A region is defined by (xmax, xmin, ymax, ymin) and the target area to
2867 be achieved for all elements in this region in the mesh specified by
2868 meshID is given by targetA. This function only performs a series of
2869 coarsenings on the elements in this region. If the area is smaller, then
2870 no refinement is done.
2871
2872 Verify correctness of the Mesh
2873 ------------------------------
2874
2875 After adaptivity operations are performed and even before adaptivity
2876 operations, it is important to first verify that we are working on a
2877 mesh that is consistent geometrically with the types of mesh that the
2878 adaptivity algorithms are designed to work on. There is a function that
2879 can be used to test various properties of a mesh, like area, quality,
2880 geometric consistency, idxl list correctness, etc.
2881
2882 ::
2883
2884   void FEM_ADAPT_TestMesh(int meshID);
2885
2886 This provides a series of tests to determine the consistency of the mesh
2887 specified by meshID.
2888
2889 These four simple steps define what needs to be used by a program that
2890 wishes to use the adaptivity features of ParFUM.
2891
2892 ParFUM developers
2893 -----------------
2894
2895 This manual is meant for ParFUM users, so developers should look at the
2896 source code and the doxygen generated documentation.
2897
2898 .. [1]
2899    The framework uses the excellent graph partitioning package Metis.