1 \documentclass[10pt]{article}
2 \usepackage{../pplmanual}
3 \input{../pplmanual}
5 \makeindex
7 \title{\charmpp\\ Iterative Finite Element Matrix (IFEM) Library\\ Manual}
8 \version{1.2}
9 \credits{
10 Initial version of \charmpp{} Finite Element Framework was developed
11 by Orion Lawlor in the spring of 2003.
12 }
14 \begin{document}
16 \maketitle
18 \section{Introduction}
20 This manual presents the Iterative Finite Element Matrix (IFEM) library,
21 a library for easily solving matrix problems derived from
22 finite-element formulations.  The library is designed to be matrix-free,
23 in that the only matrix operation required is matrix-vector product,
24 and hence the entire matrix need never be assembled
26 IFEM is built on the mesh and communication capabilities of the Charm++
27 FEM Framework, so for details on the basic runtime, problem setup, and
28 partitioning see the FEM Framework manual.
31 \subsection{Terminology}
33 A FEM program manipulates elements and nodes. An \term{element} is a portion of
34 the problem domain, also known as a cell, and is typically some simple shape
35 like a triangle, square, or hexagon in 2D; or tetrahedron or rectangular solid in 3D.
36 A \term{node} is a point in the domain, and is often the vertex of several elements.
37 Together, the elements and nodes form a \term{mesh}, which is the
38 central data structure in the FEM framework.  See the FEM manual for details.
40 % --------------- Solvers -----------------
41 \section{Solvers}
42 \label{sec:solver}
43 A IFEM \term{solver} is a subroutine that controls the search for the solution.
45 Solvers often take extra parameters, which are listed in a type called
46 in C \kw{ILSI\_Param}, which in Fortran is an array of \kw{ILSI\_PARAM} doubles.
47 You initialize these solver parameters using the subroutine \kw{ILSI\_Param\_new},
48 which takes the parameters as its only argument.  The input and output
49 parameters in an \kw{ILSI\_Param} are listed in Table~\ref{table:solverIn}
50 and Table~\ref{table:solverOut}.
52 \begin{table}[hh]
53 \begin{center}
54 \begin{tabular}{|l|l|l|}\hline
55 C Field Name & Fortran Field Offset & Use\\\hline
56 maxResidual & 1 & If nonzero, termination criteria: magnitude of residual. \\
57 maxIterations & 2 & If nonzero, termination criteria: number of iterations. \\
58 solverIn & 3-10 & Solver-specific input parameters. \\
59 \hline
60 \end{tabular}
61 \end{center}
62 \caption{\kw{ILSI\_Param} solver input parameters.}
63 \label{table:solverIn}
64 \end{table}
66 \begin{table}[hh]
67 \begin{center}
68 \begin{tabular}{|l|l|l|}\hline
69 C Field Name & Fortran Field Offset & Use\\\hline
70 residual & 11 & Magnitude of residual of final solution. \\
71 iterations & 12 & Number of iterations actually taken. \\
72 solverOut & 13-20 & Solver-specific output parameters. \\
73 \hline
74 \end{tabular}
75 \end{center}
76 \caption{\kw{ILSI\_Param} solver output parameters.}
77 \label{table:solverOut}
78 \end{table}
83 The only solver currently written using IFEM is the conjugate gradient solver.
84 This linear solver requires the matrix to be real, symmetric and positive definite.
86 Each iteration of the conjugate gradient solver requires one matrix-vector product
87 and two global dot products.  For well-conditioned problems, the solver typically
88 converges in some small multiple of the diameter of the mesh--the number of elements
89 along the largest side of the mesh.
91 You access the conjugate gradient solver via the subroutine name \kw{ILSI\_CG\_Solver}.
94 % -------------- Shared-node ----------------
95 \section{Solving Shared-Node Systems}
97 Many problems encountered in FEM analysis place the entries of the
98 known and unknown vectors at the nodes--the vertices of the domain.
99 Elements provide linear relationships between the known and unknown node values,
100 and the entire matrix expresses the combination of all these element relations.
102 For example, in a structural statics problem, we know the net force at
103 each node, $f$, and seek the displacements of each node, $u$.  Elements
104 provide the relationship, often called a stiffness matrix $K$, between
105 a nodes' displacments and its net forces:
107 $108 f=K u 109$
111 We normally label the known vector $b$ (in the example, the forces), the unknown
112 vector $x$ (in the example, the displacments), and the matrix $A$:
114 $115 b=A x 116$
119 IFEM provides two routines for solving problems of this type.  The first routine, \kw{IFEM\_Solve\_shared}, solves for the entire $x$ vector based on the known values of the $b$ vector.  The second, \kw{IFEM\_Solve\_shared\_bc}, allows certain entries in the $x$ vector to be given specific values before the problem is solved, creating values for the $b$ vector.
122 \newpage
123 \subsection{IFEM\_Solve\_shared}
125 \prototype{IFEM\_Solve\_shared}
126 \function{void IFEM\_Solve\_shared(ILSI\_Solver s,ILSI\_Param *p,
127         int fem\_mesh, int fem\_entity,int length,int width,
128         IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
129 \function{subroutine IFEM\_Solve\_shared(s,p,
130         fem\_mesh,fem\_entity,length,width,
131         A,ptr,b,x)}
132         \args{external solver subroutine :: s}
133         \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
134         \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
135         \args{external matrix-vector product subroutine :: A}
136         \args{TYPE(varies), pointer :: ptr }
137         \args{double precision, intent(in) :: b(width,length)}
138         \args{double precision, intent(inout) :: x(width,length)}
140 This routine solves the linear system $A x = b$ for the unknown vector $x$.  \uw{s} and \uw{p} give the particular linear solver to use, and are described in more detail in Section~\ref{sec:solver}.  \uw{fem\_mesh} and \uw{fem\_entity} give the FEM framework mesh (often \kw{FEM\_Mesh\_default\_read()}) and entity (often \kw{FEM\_NODE}) with which the known and unknown vectors are listed.
142 \uw{width} gives the number of degrees of freedom (entries in the vector) per node. For example, if there is one degree of freedom per node, width is one.  \uw{length} should always equal the number of FEM nodes.
144 \uw{A} is a local matrix-vector product routine you must write.  Its interface is described in Section~\ref{sec:mvp}. \uw{ptr} is a pointer passed down to \uw{A}--it is not otherwise used by the framework.
146 \uw{b} is the known vector.  \uw{x}, on input, is the initial guess for the unknown vector.  On output, \uw{x} is the final value for the unknown vector.  \uw{b} and \uw{x} should both have length * width entries.  In C, DOF $i$ of node $n$ should be indexed as $x[n*$\uw{width}$+i]$.  In Fortran, these arrays should be allocated like \uw{x(width,length)}.
148 When this routine returns, \uw{x} is the final value for the unknown vector, and the output values of the solver parameters \uw{p} will have been written.
150 \begin{alltt}
151 // C++ Example
153   int nNodes=FEM_Mesh_get_length(mesh,FEM_NODE);
154   int width=3; //A 3D problem
155   ILSI_Param solverParam;
156   struct myProblemData myData;
158   double *b=new double[nNodes*width];
159   double *x=new double[nNodes*width];
160   ... prepare solution target b and guess x ...
162   ILSI_Param_new(&solverParam);
163   solverParam.maxResidual=1.0e-4;
164   solverParam.maxIterations=500;
166   IFEM_Solve_shared(IFEM_CG_Solver,&solverParam,
167          mesh,FEM_NODE, nNodes,width,
168          myMatrixVectorProduct, &myData, b,x);
170 ! F90 Example
171   include 'ifemf.h'
172   INTEGER :: mesh, nNodes,width
173   DOUBLE PRECISION, ALLOCATABLE :: b(:,:), x(:,:)
174   DOUBLE PRECISION :: solverParam(ILSI_PARAM)
175   TYPE(myProblemData) :: myData
178   nNodes=FEM_Mesh_get_length(mesh,FEM_NODE)
179   width=3   ! A 3D problem
181   ALLOCATE(b(width,nNodes), x(width,nNodes))
182   ... prepare solution target b and guess x ..
184   ILSI_Param_new(&solverParam);
185   solverParam(1)=1.0e-4;
186   solverParam(2)=500;
188   IFEM_Solve_shared(IFEM_CG_Solver,solverParam,
189          mesh,FEM_NODE, nNodes,width,
190          myMatrixVectorProduct, myData, b,x);
192 \end{alltt}
194 \newpage
195 \subsubsection{Matrix-vector product routine}
196 \label{sec:mvp}
198 IFEM requires you to write a matrix-vector product routine that will evaluate $A x$ for various vectors $x$.  You may use any subroutine name, but it must take these arguments:
200 \prototype{IFEM\_Matrix\_product}
201 \function{void IFEM\_Matrix\_product(void *ptr,int length,int width,
202         const double *src, double *dest);}
203 \function{subroutine IFEM\_Matrix\_product(ptr,length,width,src,dest)}
204   \args{TYPE(varies), pointer  :: ptr}
205   \args{integer, intent(in)  :: length,width}
206   \args{double precision, intent(in) :: src(width,length)}
207   \args{double precision, intent(out) :: dest(width,length)}
209 The framework calls this user-written routine when it requires a matrix-vector product.  This routine should compute $dest = A \, src$, interpreting $src$ and $dest$ as vectors.  \uw{length} gives the number of nodes and \uw{width} gives the number of degrees of freedom per node, as above.
211 In writing this routine, you are responsible for choosing a representation for the matrix $A$. For many problems, there is no need to represent $A$ explicitly--instead, you simply evaluate $dest$ by looping over local elements, taking into account the values of $src$.  This example shows how to write the matrix-vector product routine for simple 1D linear elastic springs, while solving for displacement given net forces.
213 After calling this routine, the framework will handle combining the overlapping portions of these vectors across processors to arrive at a consistent global matrix-vector product.
215 \begin{alltt}
216 // C++ Example
217 #include "ifemc.h"
219 typedef struct \{
220   int nElements; //Number of local elements
221   int *conn; // Nodes adjacent to each element: 2*nElements entries
222   double k; //Uniform spring constant
223 \} myProblemData;
225 void myMatrixVectorProduct(void *ptr,int nNodes,int dofPerNode,
226           const double *src,double *dest)
227 \{
228   myProblemData *d=(myProblemData *)ptr;
229   int n,e;
230   // Zero out output force vector:
231   for (n=0;n<nNodes;n++) dest[n]=0;
232   // Add in forces from local elements
233   for (e=0;e<d->nElements;e++) {
234     int n1=d->conn[2*e+0]; // Left node
235     int n2=d->conn[2*e+1]; // Right node
236     double f=d->k * (src[n2]-src[n1]); //Force
237     dest[n1]+=f;
238     dest[n2]-=f;
239   }
240 \}
243 ! F90 Example
244   TYPE(myProblemData)
245     INTEGER :: nElements
246     INTEGER, ALLOCATABLE :: conn(2,:)
247     DOUBLE PRECISION :: k
248   END TYPE
250 SUBROUTINE myMatrixVectorProduct(d,nNodes,dofPerNode,src,dest)
251   include 'ifemf.h'
252   TYPE(myProblemData), pointer :: d
253   INTEGER :: nNodes,dofPerNode
254   DOUBLE PRECISION :: src(dofPerNode,nNodes), dest(dofPerNode,nNodes)
255   INTEGER :: e,n1,n2
256   DOUBLE PRECISION :: f
258   dest(:,:)=0.0
259   do e=1,d%nElements
260     n1=d%conn(1,e)
261     n2=d%conn(2,e)
262     f=d%k * (src(1,n2)-src(1,n1))
263     dest(1,n1)=dest(1,n1)+f
264     dest(1,n2)=dest(1,n2)+f
265   end do
266 END SUBROUTINE
267 \end{alltt}
270 \newpage
271 \subsection{IFEM\_Solve\_shared\_bc}
273 \prototype{IFEM\_Solve\_shared\_bc}
274 \function{void IFEM\_Solve\_shared\_bc(ILSI\_Solver s,ILSI\_Param *p,
275         int fem\_mesh, int fem\_entity,int length,int width,
276         int bcCount, const int *bcDOF, const double *bcValue,
277         IFEM\_Matrix\_product\_c A, void *ptr, const double *b, double *x);}
278 \function{subroutine IFEM\_Solve\_shared\_bc(s,p,
279         fem\_mesh,fem\_entity,length,width,
280         bcCount,bcDOF,bcValue,
281         A,ptr,b,x)}
282         \args{external solver subroutine :: s}
283         \args{double precision, intent(inout) :: p(ILSI\_PARAM)}
284         \args{integer, intent(in) :: fem\_mesh, fem\_entity, length,width}
285         \args{integer, intent(in) :: bcCount}
286         \args{integer, intent(in) :: bcDOF(bcCount)}
287         \args{double precision, intent(in) :: bcValue(bcCount)}
288         \args{external matrix-vector product subroutine :: A}
289         \args{TYPE(varies), pointer :: ptr }
290         \args{double precision, intent(in) :: b(width,length)}
291         \args{double precision, intent(inout) :: x(width,length)}
293 Like \kw{IFEM\_Solve\_shared}, this routine solves the linear system $A x = b$ for the unknown vector $x$.  This routine, however, adds support for boundary conditions associated with $x$. These so-called "essential" boundary conditions restrict the values of some unknowns. For example, in structural dynamics, a fixed displacment is such an essential boundary condition.
295 The only form of boundary condition currently supported is to impose a fixed value on certain unknowns, listed by their degree of freedom--that is, their entry in the unknown vector.  In general, the $i$'th DOF of node $n$ has DOF number $n*width+i$ in C and $(n-1)*width+i$ in Fortran.  The framework guarantees that, on output, for all $bcCount$ boundary conditions, $x(bcDOF(f))=bcValue(f)$.
297 For example, if $width$ is 3 in a 3d problem, we would set node $ny$'s y coordinate to 4.6 and node $nz$'s z coordinate to 7.3 like this:
299 \begin{alltt}
300 // C++ Example
301   int bcCount=2;
302   int bcDOF[bcCount];
303   double bcValue[bcCount];
304   // Fix node ny's y coordinate
305   bcDOF=ny*width+1; // y is coordinate 1
306   bcValue=4.6;
307   // Fix node nz's z coordinate
308   bcDOF=nz*width+2; // z is coordinate 2
309   bcValue=2.0;
311 ! F90 Example
312 // C++ Example
313   integer :: bcCount=2;
314   integer :: bcDOF(bcCount);
315   double precision :: bcValue(bcCount);
316   // Fix node ny's y coordinate
317   bcDOF(1)=(ny-1)*width+2; // y is coordinate 2
318   bcValue(1)=4.6;
319   // Fix node nz's z coordinate
320   bcDOF(2)=(nz-1)*width+3; // z is coordinate 3
321   bcValue(2)=2.0;
322 \end{alltt}
326 Mathematically, what is happening is we are splitting the partially unknown vector $x$ into a completely unknown portion $y$ and a known part $f$:
327 $A x = b$
328 $A (y + f) = b$
329 $A y = b - A f$
331 We can then define a new right hand side vector $c=b-A f$ and solve the new linear system $A y=c$ normally.  Rather than renumbering, we do this by zeroing out the known portion of $x$ to make $y$.  The creation of the new linear system, and the substitution back to solve the original system are all done inside this subroutine.
333 One important missing feature is the ability to specify general linear constraints on the unknowns, rather than imposing specific values.
336 \section{Index}
337 \input{index}
338 \end{document}