Synopsis:

• &t peforms the summation convention multiplication of Cartesian tensors. It is a binary operator that looks for indices in the two operands, and whenever it finds two indices that are equal it performs a sum from 1 to 3 over that index.
• &t can be combined with the nab operator and the predefined tensors delta[i,j] and epsilon[i,j,k] to produce cross products, gradients, curls etc.
• &t follows the usual arithemetic rules, parentheses can be used etc.
• The user must take care that the indices over which summation should be performed have not been assigned values.
• There is little or no syntax check. The expressions a[i] &t b[i,k] &t c[i,j] and a[i] &t ( s + b[i]) are meaningless, but &t will not give any explicit warning.

Examples:

> read(femLego);
> a[i] &t b[i];
       a[1] b[1] + a[2] b[2] + a[3] b[3]

> s &t a[k];
                                     s a[k]

> r &t s;
                                      r s

> u[i](x,y,z) &t nab(u[j](x,y,z))[i];

                  /  d               \                 /  d               \
    u[1](x, y, z) |---- u[j](x, y, z)| + u[2](x, y, z) |---- u[j](x, y, z)|
                  \ dx               /                 \ dy               /

                         /  d               \
         + u[3](x, y, z) |---- u[j](x, y, z)|
                         \ dz               /

> epsilon[1,j,k] &t a[j] &t b[k];

                             a[2] b[3] - a[3] b[2]

> for i from 1 to 3 do   epsilon[i,j,k] &t a[j] &t b[k]   od;

                             a[2] b[3] - a[3] b[2]

                            - a[1] b[3] + a[3] b[1]

                             a[1] b[2] - a[2] b[1]

nab(f)[i]

Parameters:

f - an expression or function depending on x,y,z.

• nab(f)[i] returns the i-th component of the gradient of the argument f. f is assumed to depend on x,y and/or z.

Examples:

> nab(u(x,y,z))[i];

                               nab(u(x, y, z))[i]

> nab(phi(x,y,z))[i] &t nab(u(x,y,z))[i] ;

 /  d              \ /  d            \   /  d              \ /  d            \
 |---- phi(x, y, z)| |---- u(x, y, z)| + |---- phi(x, y, z)| |---- u(x, y, z)|
 \ dx              / \ dx            /   \ dy              / \ dy            /

        /  d              \ /  d            \
      + |---- phi(x, y, z)| |---- u(x, y, z)|
        \ dz              / \ dz            /

> delta[i,j] &t nab(   nab(u(x,y,z))[i] )[j];

          /   2            \   /   2            \   /   2            \
          |  d             |   |  d             |   |  d             |
          |----- u(x, y, z)| + |----- u(x, y, z)| + |----- u(x, y, z)|
          |   2            |   |   2            |   |   2            |
          \ dx             /   \ dy             /   \ dz             /

> delta[i,j] &t nab(   nab(u(x,y))[i] )[j];

                       /   2         \   /   2         \
                       |  d          |   |  d          |
                       |----- u(x, y)| + |----- u(x, y)|
                       |   2         |   |   2         |
                       \ dx          /   \ dy          /

Synopsis:

• The unit Cartesian tensor (delta[i,j] = 1 if i=j, 0 otherwise). To be used with &t and nab(.)[i] to form general Cartesian tensor expressions.

Examples:

> delta[1,1];
                                       1
> delta[1,2];
                                       0
> A[i,j] &t delta[j,k] &t b[k,n];
              A[i, 1] b[1, n] + A[i, 2] b[2, n] + A[i, 3] b[3, n]

Synopsis:

• The Cartesian permutation tensor. To be used with &t and nab(.)[i] to form general Cartesian tensor expressions.

Examples:

> epsilon[1,1,1];
                                       0
> epsilon[1,2,3];
                                       1
> epsilon[1,3,2];
                                       -1
> epsilon[1,j,k] &t a[j] &t b[k];
                             a[2] b[3] - a[3] b[2]
> for i from 1 to 3 do  epsilon[i,j,k] &t a[j] &t b[k]  od; i:='i':
                             a[2] b[3] - a[3] b[2]
                            - a[1] b[3] + a[3] b[1]
                             a[1] b[2] - a[2] b[1]

partInt(expr,var):

Parameters

expr - The expression that is to be partially integrated.

var - variable or list of variables appearing inside ElementInt() in expr.

• The expression in expr is integrated partially once, so that the order of diffferentiation of the variable(s) in var is reduced by one. expr is expected to contain ElementInt() integrals, only such are affected.
• var can be either a variable appearing in expression or a list of variables. If var is a list, one partial integration is performed for each variable, so that partInt(expr,[u(x),v(x)]) is equivalent to partInt(partInt(expr,u(x)),v(x)).
• The partial integrations produce boundary integrals, BoundaryInt(), and components of normal vectors, called BdNx, BdNy etc. The normal vector is expected to be outwardly directed and normal to the boundary.
• When expressions generated by partInt are passed to the code generation commands, such as mkFem, the normal vector components BdNx, BdNy etc., are replaced by their appropriate definitions which should appear in the *_bsf files specifying the elements.

Examples:

> partInt(ElementInt(f(x)*diff(g(x),x)),g(x));

                 /d      \
-ElementInt(g(x) |-- f(x)|) + BoundaryInt(BdNx g(x) f(x))
                 \dx     /

> partInt( 0 = ElementInt(f(x)*diff(g(x),x)) ,g(x));

                         /d      \
    0 = -ElementInt(g(x) |-- f(x)|) + BoundaryInt(BdNx g(x) f(x))
                         \dx     /

> Lapl:= nab(nab(u(x,y))[i])[j] &t delta[i,j];

                        / 2         \   / 2         \
                        |d          |   |d          |
                Lapl := |--- u(x, y)| + |--- u(x, y)|
                        |  2        |   |  2        |
                        \dx         /   \dy         /

> partInt( ElementInt(test(x,y)*Lapl) ,u(x,y));

            /d         \ /d            \
-ElementInt(|-- u(x, y)| |-- test(x, y)|)
            \dx        / \dx           /

                        /d         \
     + BoundaryInt(BdNx |-- u(x, y)| test(x, y))
                        \dx        /

                  /d         \ /d            \
     - ElementInt(|-- u(x, y)| |-- test(x, y)|)
                  \dy        / \dy           /

                        /d         \
     + BoundaryInt(BdNy |-- u(x, y)| test(x, y))
                        \dy        /

BCsubs(bclist,expr):

Parameters

bclist - A list of equalities specifying boundary conditions. The left hand sides must be either variables, or derivatives of variables, appearing inside BoundaryInt() in expr.

expr - An expression containing BoundaryInt()s, typically as a result of using partInt.

• In deriving a variational form of a PDE, Neumann boundary conditions are typically introduced by doing the appropriate substitutions in the boundary integrals that appear when the expression is integrated partially. BCsubs is a tool for doing this.
• BCsubs works in a way similar to subs. Occurencies of the left hand sides of the equalitites in BClist are replaced by the right hand sides.
• BCsubs will only affect arguments of BoundaryInt. All other terms are returned unchanged.
• If the left hand side of an entry in BClist is an undifferentiated variable, only undifferentiated occurencies of this variable are replaced. Similarly a first derivative will only replace first derivatives, etc. This is the main difference between BCsubs and subs.

Examples:

> femform:=0=partInt(ElementInt(test(x)*diff(v(x),x,x)),v(x));

                           /d      \ /d         \
femform := 0 = -ElementInt(|-- v(x)| |-- test(x)|)
                           \dx     / \dx        /

                        /d      \
     + BoundaryInt(BdNx |-- v(x)| test(x))
                        \dx     /

> BCsubs([diff(v(x),x)=q],femform);

                /d      \ /d         \
0 = -ElementInt(|-- v(x)| |-- test(x)|) + BoundaryInt(BdNx q test(x))
                \dx     / \dx        /

variation(expr,var,dvar)

Parameters

expr - An expression containing the variables in var.

var - A variable, or a list of variables, to perform the variation with respect to.

dvar - the names of the variations of the variables in var. Must have same length as var

• Performs the variation of expression by replacing var by var + eps*dvar, differentiating with repect to eps, and substituting eps=0.
• ElementInt and BoundaryInt commute with this differentiation, so that variations of integral forms expressed in ElementInt and BoundaryInt can be computed.
• if var is a list, the chain rule is applied so that variation(expr,[u(x),v(x)],[du(x),dv(x)]) is equivalent to
  variation(expr,u(x),du(x)) + variation(expr,v(x),dv(x))
• expr can be a list, equality or a plain expression.

Examples:

> J:=ElementInt(u(x)^2);
                                            2
                        J := ElementInt(u(x) )

> variation( J , u(x),du(x));

                       ElementInt(2 u(x) du(x))


> F:=ElementInt(u(x,y)^2 + v(x,y)^2);
                                        2          2
                 F := ElementInt(u(x, y)  + v(x, y) )

> variation(F,[u(x,y),v(x,y)],[du(x,y),dv(x,y)]);

   ElementInt(2 u(x, y) du(x, y)) + ElementInt(2 v(x, y) dv(x, y))


> BClist:=[u(x)=0,diff(u(x),x)=q];

                                       d
                  BClist := [u(x) = 0, -- u(x) = q]
                                       dx

> variation(BClist,u(x),du(x));

                                  d
                      [du(x) = 0, -- du(x) = 0]
                                  dx

test_subs(test_list,testfcn,expr)

Parameters

test_list - a list of test functions

test_fcn - a name of a variable to be used as a test function

expr - An expression.

• When deriving FEM formulations for systems of PDEs, it may be convenient to multiply each equation by a different test function, and add them together to one large integral form. To recover the individual equations in a form suitable for input to for example mkFem, all test functions except one are replaced by 0. test_subs provides a convenient shorthand for this.
• test_list is a list of the names of the different formal test functions.
• testfcn is a name of a test function to be used with for example mkFem.
• The output is a list of equations in a form suitable for input to mkFem.

Examples:

> ueq:=0=diff(u(x),x,x);

                                     2
                                    d
                         ueq := 0 = --- u(x)
                                      2
                                    dx

> veq:=u(x)=diff(v(x),x,x);

                                       2
                                      d
                        veq := u(x) = --- v(x)
                                        2
                                      dx

> varForm:=ElementInt(ut(x) * (lhs(ueq)-rhs(ueq)))+
           ElementInt(vt(x) * (lhs(veq)-rhs(veq))) =0:

> pivarForm:=partInt(varForm,[u(x),v(x)]):
> BCpivarForm:=BCsubs([diff(u(x),x)=q,diff(v(x),x)=0],pivarForm);

                          /d      \ /d       \
BCpivarForm := ElementInt(|-- u(x)| |-- ut(x)|)
                          \dx     / \dx      /

     - BoundaryInt(BdNx q ut(x)) + ElementInt(vt(x) u(x))

                  /d      \ /d       \
     + ElementInt(|-- v(x)| |-- vt(x)|) = 0
                  \dx     / \dx      /

> eqns:=test_subs([ut(x),vt(x)],test(x),BCpivarForm):
> eqns[1];eqns[2];

            /d      \ /d         \
 ElementInt(|-- u(x)| |-- test(x)|) - BoundaryInt(BdNx q test(x)) = 0
            \dx     / \dx        /


                                        /d      \ /d         \
  ElementInt(test(x) u(x)) + ElementInt(|-- v(x)| |-- test(x)|) = 0
                                        \dx     / \dx        /

getFemLabInput(params_list):

Parameters

params_list - A list of scalar input parameters.

• Generates the code to read meshes generated by FEMLAB 1.1. When the solver is run it expects to find a file called area.scratch containing the mesh.
• getFemLabInput also generates code to read the scalar parameters that are specified in params_list from standard input. getFemLabInput writes a file called indata_sample in the work directory, which can be used as a template for a file specifying all input parameters. After editing indata_sample , the parameter values are conveniently passed to the solver by the command a.out < indata_sample.

Examples:

> # get files to read input from femlab;
> getFemLabInput(params_list)
rdinpt.f
rdpara.f
indata_sample

getFemLabInputP2P1(params_list):

Parameters

params_list - A list of scalar input parameters.

• This modifies the mesh information for 2D triangle meshes from FEMLAB 1.1 so that it can be used with the basefunction file 2DP2P1_gsq_bsf. The function is the same as getFemLabInput otherwise.
• Generates the code to read meshes generated by FEMLAB 1.1. When the solver is run it expects to find a file called area.scratch containing the mesh.
• getFemLabInput also generates code to read the scalar parameters that are specified in params_list from standard input. getFemLabInput writes a file called indata_sample in the work directory, which can be used as a template for a file specifying all input parameters. After editing indata_sample , the parameter values are conveniently passed to the solver by the command a.out < indata_sample.

Examples:

> # get files to read input from femlab;
> getFemLabInput(params_list)
rdinpt.f
rdpara.f
indata_sample

FemLabPlot(plot_list,menu_list, unknown_list);

Parameters

plot_list - A list of items to plot. If an item is a list of two variables, it will be plotted as a vector field by FEMLAB.

menu_list - A list of strings to appear in the FEMLAB menus. There should be one entry for each entry in plot_list.

unknown_list - The list of unknowns as discussed above.

• Creates code to save computed results in a form readable by FEMLAB. Output will be saved in files named outN.d, where N is an increasing integer value.
• Creates the file .femlab_modules which informs FEMLAB about what the variables are, menu appearances, etc.
• An entry in 'plot_list' that is an unknown from 'unknown_list' will cause that variable to be plotted as a scalar field. If an entry in 'plot_list' is a list of two variables from 'unknown_list', this pair will be plotted as a vector field.
• FemLabPlot will work also with the mixed element 2DP2P1_gsq_bsf, by storing only nodal values at mesh junctions and discarding nodal values at triangle edges.

Examples:

> # plot u(x,y), v(x,y) as two scalar fields;
> FemLabPlot([u(x,y) ,v(x,y)],[u,v],unknown_list);
.femlab_modules
wroutp
sv2fml

> # plot u(x,y),v(x,y) as a vector field and p(x,y) as a scalar;
> FemLabPlot([[u(x,y),v(x,y)], p(x,y)] [uv_vector,pressure],unknown_list);
.femlab_modules
wroutp
sv2fml

get3DP1_input(params_list);

Parameters:

params_list - A list of scalar input parameters.

• Generates the code to read 3D tetrahedral meshes, as created for instance by mk3DP1_cube and mk3DP1_bcond. When the solver is run it expects to find a file called mesh.scratch containing the mesh, and a file called bcond.scratch containing boundary conditions.
• get3DP1_input also generates code to read the scalar parameters that are specified in params_list from standard input. get3DP1_input creates a file called indata_sample in the work directory, which can be used as a template for a file specifying all input parameters. After editing indata_sample , the parameter values are conveniently passed to the solver by the command a.out < indata_sample.
• This is included as a demonstration of the 3D element. A coupling to a serious 3D mesh generator is still missing.

Examples:

> # set up to read input data: copies rdinpt.f to ./source;
> get3DP1_input(params_list);
rdinpt.f
rdpara.f
indata_sample

mk3DP1_cube(box,npts);

Parameters:

box - a list of the form [xlo, xhi, ylo, yhi, zlo, zhi], where xlo - zhi are numeric values that specify the range of the three coordinates x,y,z.

npts - a list with three integer entries, specifying the number of nodes in the x,y and z directions.

• Make a uniform 3D grid in a rectangular box. The mesh is saved in the file mesh.scratch.
• After completion the number of nodes and elements are echoed.
• This is included as a demonstration of the 3D element. A coupling to a serious 3D mesh generator is still missing.

Examples:

> # make a 10x10x10 mesh in the unit cube:
> # box:=[xlo, xhi, ylo, yhi, zlo, zhi];
> box:=[0.,1., 0.,1., 0.,1.];
> npts:=[10,10,10];
> mk3DP1_cube(box,npts);
[1331,6000]

mk3DP1_bcond(bc_list);

Parameters:

bc_list - a list of the form [u_bc, v_bc, ..], with one entry for each dependent variable in the problem. Each of the u_bc have the form [[c1,r1], [c2,r2], [c3,r3], [c4,r4], [c5,r5], [c6,r6]], There is one pair [c, r] for each face of the block. The first member can be either D or N: c=N will give Neumann boundary conditions, c=D: Dirichlet. The second member r should be a numerical value which is made available to the solver as qBCval. The faces are numbered acccording to 1: y = ylo, 2: x=xhi, 3: y=yhi, 4: x=xlo, 5: z=zlo, 6: z=zhi.

• Supply boundary condition information for uniform 3D grids generated by mk3DP1_cube. The information is saved in the file bcond.scratch.
• This is included as a demonstration of the 3D element. A coupling to a serious 3D mesh generator is still missing.

> # make boundary conditions file (bcond.scratch);
> # supply a pair [c, r] for each face of the block:
> # c=N: Neumann BC. The numerical value r is available as qBCval.
> # c=D: Dirichlet BC. The numerical value r is available as qBCval.
> # The faces are numbered acccording to
# 1: y = ylo, 2: x=xhi, 3: y=yhi, 4: x=xlo, 5: z=zlo, 6: z=zhi;
> u_bc:=[[N,0.],[N,0.],[N,0.],[N,0.],[D,0.],[N,1.]];
> v_bc:=[[D,0.],[D,0.],[D,0.],[D,0.],[D,0.],[N,1.]];
> mk3DP1_bcond( [u_bc, v_bc] );

AVSPlot_3DP1(unknown_list);

Parameters

unknown_list - The list of unknowns as discussed above.

• Set up code to save results in a form readable by AVS. Output will be saved in files named outN.inp, where N is an increasing integer value.
• Should only be used with 3D tetrahedral elements, as 3DP1_gsq_bsf.
• This is included as a demonstration of the 3D element. However, when a coupling to a serious 3D mesh generator is made, AVSPlot_3DP1 should work without change.

> # set up to save results to plot with AVS:
> AVSPlot_3DP1(unknown_list);
wroutp

mkDirBC(DBCeq, unknown_list, old_unknown_list, params_list, bsf_file)

Parameters:

DBCeq - A list of equations specifying Dirichlet boundary conditions. Each entry in the list is an equation of the form u(x,y)=expr, where u(x,y) is a variable in unknown_list. For those variables in unknown_list, which do not appear in DBCeq, a default Dirichlet boundary condition is used. This is u(x,y) = qBCval, where qBCval is a predefined value which is specified with the mesh.
If one of the entries in DBCeq is SinglePoint(var), where var is one of the variables in unknown_list, a Dirichlet condition is applied at one single node, prescribing the variable to be zero at that node.

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.

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

• Generates the code needed to specify Dirichlet boundary conditions.
• If the string 'qBCval' is encountered in DBCeq, it is replaced by a predefined value specified with the mesh. This gives a way of separating parts of the boundary with different Dirichlet conditions.
• If one of the entries in DBCeq is SinglePoint(var), the value of that variable will be set to zero at one single node. This is useful for equations with only Neumann conditons, for which the solution is only determined up to an additive constant. Using SinglePoint(var) this additive constant is fixed (arbitrarily), and possible numerical problems are avoided.

Examples:

> # set up lists:
> eq_list:= [ueq,veq,peq]:
> unknown_list:= [u(x,y),v(x,y),p(x,y)]:
> old_unknown_list:= [uo(x,y),vo(x,y),po(x,y)]:
> params_list:= [Diffusivity,bigQ,BCpara]:
> # Dirichlet BCs;
> # The call below creates the following boundary conditions:
# u(x,y)=x*BCpara on Dirichlet parts of the boundary, ;
# v(x,y)=qBCval on Dirichlet parts of the boundary;
# p(x,y)=0 at one single node;
> DBCueq:= u(x,y)=x*BCpara:
> mkDirBC([DBCueq,SinglePoint(p(x,y))],unknown_list,
old_unknown_list,params_list,`2DP1_gsq_bsf`):

mkDirBC2(DBCeq, unknown_list, old_unknown_list, params_list,unkntestlist,testbsflist)

Parameters:

DBCeq - A list of equations specifying Dirichlet boundary conditions. Each entry in the list is an equation of the form u(x,y)=expr, where u(x,y) is a variable in unknown_list. For those variables in unknown_list, which do not appear in DBCeq, a default Dirichlet boundary condition is used. This is u(x,y) = qBCval, where qBCval is a predefined value which is specified with the mesh.
If one of the entries in DBCeq is SinglePoint(var), where var is one of the variables in unknown_list, a Dirichlet condition is applied at one single node, prescribing the variable to be zero at that node.

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.

unkntestlist - A list of testfunction names that couple a certain testfunction with a certain unknown. The number of entries in unkntestlist must be the same as in unknown_list and old_unknown_list. The usage is the same as for the corresponding argument to mkFem2.

testbsflist - A list containing three entries, the first two being names of testfunctions, the last the name of a file that contins the definitions of the basefunctions. The two basefunction names are identified with the actual basefunctions via ordering and the name of this file, i.e., with testbsflist:=[v(x,y),q(x,y),`2DP2P1_gsq_bsf`], v(x,y) is P2 and q(x,y) is P1. The usage is the same as for the corresponding argument to mkFem2.

• Generates the code needed to specify Dirichlet boundary conditions for use with mixed formulations with two different basefunctions. Except for the last two arguments the function is the same as mkDirBC.
• If the string 'qBCval' is encountered in DBCeq, it is replaced by a predefined value specified with the mesh. This gives a way of separating parts of the boundary with different Dirichlet conditions.
• If one of the entries in DBCeq is SinglePoint(var), the value of that variable will be set to zero at one single node. This is useful for equations with only Neumann conditons, for which the solution is only determined up to an additive constant. Using SinglePoint(var) this additive constant is fixed (arbitrarily), and possible numerical problems are avoided.

Examples:

> # set up lists:
> eq_list:= [u1eq,u2eq,peq]:
> unknown_list:= [u1(x,y),u2(x,y),p(x,y)]:
> old_unknown_list:= [uo1(x,y),uo2(x,y),po(x,y)]:
> params_list:=[Diffusivity,bigQ,BCpara]:
> unkntestlist:=[v(x,y),v(x,y),q(x,y)]:
> testbsflist:=[v(x,y),q(x,y),`2DP2P1_gsq_bsf`]:
> # The unknowns u1(x,y),u2(x,y) are associated with
> # the test function v(x,y), and p(x,y) with q(x,y).
> # v(x,y) represents P2 basefunctions, and q(x,y) represents P1.
>
> # Dirichlet BCs;
> # The call below creates the following boundary conditions:
# u(x,y)=x*BCpara on Dirichlet parts of the boundary, ;
# v(x,y)=qBCval on Dirichlet parts of the boundary;
# p(x,y)=0 at one single node;
> DBCueq:= u(x,y)=x*BCpara:
> mkDirBC([DBCueq,SinglePoint(p(x,y))],unknown_list,
old_unknown_list,params_list,`2DP1_gsq_bsf`):

mkIC(ICeq, unknown_list, params_list, bsf_file);

Parameters:

ICeq - A list of equations specifying initial conditions. Each entry in the list is an equation of the form u(x,y)=expr, where u(x,y) is a variable in unknown_list. Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

unknown_list - The list of unknowns as discussed above.

params_list - A list of scalar input parameters.

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

• Generates the code needed to set initial conditions. During execution the names of generated subroutines are echoed. The equations in ICeq are projected onto the test functions and the initial conditions are solved for using conjugate gradients.
• Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

Examples:

> # make initial conditions. ic_eq is on the form 'unknown = expr':
> # variables that do not enter in lhs of an ic_eq are given the
# intial value 0:
> mkIC([ v(x,y)=x*x+y*y ], unknown_list, params_list ,
`2DP1_gsq_bsf`);
ic_init
icmkre
icamdp
icadr1
icadm1
icadr2
icadm2

mkICcopy(ICeq, unknown_list, params_list, bsf_file);

Parameters:

ICeq - A list of equations specifying initial conditions. Each entry in the list is an equation of the form u(x,y)=expr, where u(x,y) is a variable in unknown_list. Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

unknown_list - The list of unknowns as discussed above.

params_list - A list of scalar input parameters.

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

• Generates the code needed to set initial conditions. The nodal values of the unknowns are assigned the values of the right hand sides of the equations in ICeq, evaluated at the corresponding points.
• Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

Examples:

mkICcopy2(ICeq, unknown_list, params_list, unkntestlist,testbsflist);

Parameters:

ICeq - A list of equations specifying initial conditions. Each entry in the list is an equation of the form u(x,y)=expr, where u(x,y) is a variable in unknown_list. Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

unknown_list - The list of unknowns as discussed above.

params_list - A list of scalar input parameters.

unkntestlist - A list of testfunction names that couple a certain testfunction with a certain unknown. The number of entries in unkntestlist must be the same as in unknown_list and old_unknown_list. The usage is the same as for the corresponding argument to mkFem2.

testbsflist - A list containing three entries, the first two being names of testfunctions, the last the name of a file that contins the definitions of the basefunctions. The two basefunction names are identified with the actual basefunctions via ordering and the name of this file, i.e., with testbsflist:=[v(x,y),q(x,y),`2DP2P1_gsq_bsf`], v(x,y) is P2 and q(x,y) is P1. The usage is the same as for the corresponding argument to mkFem2.

• Generates the code needed to set initial conditions. The nodal values of the unknowns are assigned the values of the right hand sides of the equations in ICeq, evaluated at the corresponding points.
• Those variables in unknown_list which do not appear on the left hand side of an equation in ICeq, are set to zero initially.

Examples:

> # make initial conditions. ic_eq is on the form 'unknown = expr':
> # variables that do not enter in lhs of an ic_eq are given the
# intial value 0:
> mkICcopy2([ v(x,y)=x*x+y*y ], unknown_list, params_list, unkntestlist, testbsflist);

ic_init

mkSolve(solve_list,eq_list)

Parameters:

solve_list - List of keywords to denote different linear systems solver. One keyword for each equation in eq_list.

eq_list - List of equations to be solved.

• Generates the calls to different linear systems solvers. The equations are (so far) solved in a split time stepping. Thus each equation is considered separately, and solved for one of the variables. The first equation in eq_list is solved using the solver specified by the first entry in solve_list, etc.
• Presently, the following solver keywords are supported:

iccg
Calls the Incomplete Cholesky preconditioned Conjugate Gradient solver in SLAP. It requires the linear system to be symmetric and positive definite.

gmres
Calls the GMRES iterative solver in SLAP. This allows the system to be nonsymmetric but is considerably slower than iccg for symmetric systems.

diag
This assumes the matrix to be diagonal, and solves the system directly. Typically this solver would be used in explicit time stepping schemes together with the Lump modifer to mkFem.

> # set up different solvers to use with different equations:
> solve_list:=[iccg,iccg]:
> mkSolve(solve_list,eq_list);
solve

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.
• 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 members of the equations in eq_list must be sums of the functions ElementInt and BoundaryInt.
• 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.
• 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 three values that should be entered in include/size.h. These depend on the finite element used.
• 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.

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.

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

mkFem2(eq_list, unknown_list, old_unknown_list, params_list, unkntestlist, testbsflist);

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.

unkntestlist - A list of testfunction names that couple a certain testfunction with a certain unknown. The number of entries in unkntestlist must be the same as in unknown_list and old_unknown_list.

testbsflist - A list containing three entries, the first two being names of testfunctions, the last the name of a file that contins the definitions of the basefunctions. The two basefunction names are identified with the actual basefunctions via ordering and the name of this file, i.e., with testbsflist:=[v(x,y),q(x,y),`2DP2P1_gsq_bsf`], v(x,y) is P2 and q(x,y) is P1.

• 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.
• 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 members of the equations in eq_list must be sums of the functions ElementInt and BoundaryInt.
• 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.
• 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 mkFem2 prints three values that should be entered in include/size.h. These depend on the finite element used. In particular the normal use of mkFem2 with a mixed formualtion requires that noofmx = 2 in include/size.h.
• 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.

ElementIntHiQ(expr)
This can be used to invoke an alternate integral evaluation, replacing ElementInt(), for selected terms, giving for instance a more accurate Gauss quadrature fomula. It relies on finding a definition of the procedure myElIntHiQ in the basefunction file.

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.

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).

Lump2(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