mkFem(eq_list, unknown_list, old_unknown_list, params_list, testfcn, bsf_file);

Parameters:

eq_list - A list of equations defining a PDE problem in variational form.

unknown_list - The list of unknowns as discussed above.

old_unknown_list - The list of unknowns at the previous timestep, as discussed above.

params_list - A list of scalar input parameters.

testfcn - The name of the test function used in the equations in eq_list.

bsf_file - A file containing a Maple description of a particular finite element.

• mkFem generates subroutines that compute residuals and matrix elements for the equations specified in eq_list. These are assumed to represent a variational formulation of a system of PDEs.
• The members of the equations in eq_list must be sums of the functions ElementInt() and BoundaryInt(), or the entire equation be enclosed in a CopyEq().
• If defined in the basefunction file (as myElIntHiQ), an additional integration operator, ElementIntHiQ() may be used. This can be used to invoke an alternate integral evaluation, for instance a more accurate Gauss quadrature fomula, for selected terms.
• To allow variables to be defined by simple algebraic expressions, an equation in eq_list may be written as CopyEq( u(x,y) = v(x,y)^2 ). In this example, the nodal values for variable u(x,y) will be assigned to the square of the corresponding nodal values of v(x,y). Derivatives are meaningless inside a CopyEq( ). The solver specified for such an equation must be set as 'copy' in mkSolve or mkSolveAdapt.
• It might be important that the algebraic system that is generated is positive definite. This will be the case if, for a well posed diffusion equation, the time derivative is on the left hand side.
• The solver will ask for numerical values of the parameters listed in params_list. The strings in params_list will thus be replaced by those values at runtime.
• The strings 'dt' and 'time' will hold the timestep dt and the actual time in the simulation. The solver will ask for dt when it is started, and 'time' will be incremented by dt once in each loop. Explicit time dependencies in equations or boundary conditions can thus be included simply by writing expressions as functions of 'time'.
• All function names that are not resolved by Maple will be passed on to the fortran subroutines. Thus it makes perfect sense to write for instance Sin(u(x,y)) inside an ElementInt() call. This will be computed as the sine of local values of the variable u(x,y), as it should. Any other function name that is defined in fortran (or supplied by the user) can be used in a similar fashion.
• When the generation is complete mkFem prints a few integer values that should be entered in include/size.h. These depend on the finite element used and the structure of the equations.
• mkFem prints 'changing jacobians list:' and a list of logical values after generation is complete: mkFem determines whether the system matrix in each equation changes during the solution, or if the first assembled matrix can be used subsequently. The list of logical values are one for each equation, in the order in eq_list, with 'false' meaning that the matrix is not changed, and 'true' that it has to be reassembled in each timestep. This should all work automatically and is provided as a check only.
• mkFem looks for function names in eq_list, which control how the finite element approximations are made. The functions that are understood presently are:

ElementInt(expr)
Generates local residuals by integration over an element. The equations in eq_list must consist of terms enclosed in ElementInt or BoundaryInt calls, or the entire equation be enclosed in a CopyEq().

BoundaryInt(expr)
The BoundaryInt call is used to evaluate the boundary integrals that appear in natural boundary conditions in variational formulations. BoundaryInt generates local residuals by integration over element sides. These are assembled for sides located on parts of the boundary that have 'Neumann' boundary conditions. If the string 'qBCval' is encountered it will be replaced by a number read with the mesh when the solver is run. The equations in eq_list must consist of terms enclosed in ElementInt or BoundaryInt calls, or the entire equation be enclosed in a CopyEq().

CopyEq(variable=expression)
'variable' is one of the unknowns in unknown_list.'expression' must not contain derivatives. The nodal values of 'variable' will be assigned to 'expression' evaluated at that node. The solver specified in mkSolve or mkSolveAdapt for the corresponding equation must be 'copy'.

NoExpand(expr)
This is meaningful only when used in the argument of an ElementInt or BoundaryInt function. It alters the way the symbolic calculations are done, but should not affect the numerical results. It can be used to prevent Maple from expanding nonlinear expressions of variables into huge, awkward expressions. For example: ElementInt(test(x,y)*(NoExpand(u(x,y)) )^3) could produce less code than ElementInt(test(x,y)*(u(x,y))^3).

Lump(expr)
This is an implementation of the familiar FEM-trick of 'lumping the mass matrix'. Typically this is useful in explicit schemes when the unknown u(x,y) is to be solved from something like ElementInt(test(x,y)*u(x,y))=RHS. 'Lumping the mass matrix' means to replace the system matrix on the left hand side by a diagonal matrix ('lumping' the nodes onto the diagonal). So if instead we write ElementInt( Lump(test(x,y)*u(x,y)) ) = RHS, this equation can be solved faster by using the diagonal solver (diag in mkSolve). Lump() is meaningful only when used in the argument of an ElementInt or BoundaryInt function.

PointWise(expr)
This is an implementation of a trick sometimes used in FEM, for cheaper evaluation of nonlinear expressions. It will apply a nonlinear expression to the nodal values directly, instead of first expanding the unknown in base functions, i.e.: u(x,y)^3 is expanded in base functions as
sum_i( u_i^3 * test_i) instead of ( sum_i( ui * test_i))^3. PointWise() is meaningful only when used in the argument of an ElementInt or BoundaryInt function.

Examples:

> # create residual computations etc:
> size_parameters:=mkFem(eq_list,unknown_list,old_unknown_list,params_list,test(x,y), `2DP1_gsq_bsf`);
mkresi
amedsp
addr1
addm1
addr2
addm2
grbcrs
bdedsp
adbr1
adbm1
adbr2
adbm2
addini


check include/size.h:
nvar: 2
nobf: 3
bnobf: 2

changing jacobians list:
false false