Suggestions for final projects

Various options

• Suggest yourself, can also be related to Cooper's etc. stuff
• Implement some part of a FEM-solver, include some interesting applications, study errors. Partly with PDEtool/Femlab
• Study time-dependent problems (parabolic, hyperbolic). Perhaps less implementation and more applications
• Combinations of (a subset of) the above items

Why FEM-solver - isn't it enough to have toolbox

• To learn to understand the method implement some parts of it
• Add new features. Toolbox/Femlab contains only linear basis functions, which can be regarded as a major limitation.
• You might think of other possibilities as well like post-proseesing tools, tools that help understand the various FEM ingredients (FEM-learning/teaching tools) etc.

Fem-solver for elliptic problems

Problem 3.16 pp. 175-176 in Celia-Grey (These 2 pages are not needed, all is here): Write a Matlab-code for Poisson's equation given the boundary conditions of Fig. 3.15 (handouts) and f(x,y)=0. Create the triangular mesh with PDEtool.

The flow of the program

• Define the domain and create the triangularization in pdetool .
• Define the equation somehow (you) .
• Assemble the linear system including the boundary conditions. Result should be a linear system of type Au=b.
• Solve the system using matlab sparse commands (simply u=A\b or use pcg).
• To visualize the results create a function that evaluates the solution on given canonical coordinates within a given triangle (=index in matrix t). (We'll see if we can do the rest)

Set up of the differential equation

div( c grad u) + d u = f

Here c,d ja f are matlab functions (R² -> R), of type

function f=foo(x,y)

Give the names of the functions as parameters to the assemby part. Use the feval function, to evaluate functions.

Elements

For the start, we recommed you to start with linear and quadratic elements. See the handouts.

Rectangular elements In spite of the above recommendation we will give some hints on working with rectangular elements. In this case you can't use the toolbox-data structures directly. On the other hand the data structures are simpler. Also, you might consider using the "poi"- "rectangular triangulation" functions.

Boundaries and boundary conditions

Fix (initial) sections of the boundary name them. Give them a parametrization (approximatively the curve length). Pdetool updates the E-matrix accordingly automatically. It even evaluates this right if given to initmesh.m

Write a separate file for evaluation of the boundary-conditions, pass the relative position on that edge as a parameter. (From this the corresponding parameter value can be computed.) The values of the function $u$ and the orthogonal derivative $u_n$ has to be passed. This way you deserve the maximal freedom of writing different boundary conditions. Also many of the difficulties of the parameter passing are circumvented.

See also pdedemo3.m demo.

Assembling the system

It is important to note that there are three different types of basis functions, classified on their continuity properties over the triangles:

• Elements associated with nodes
• Elements associated with sides of triangles
• Elements associated with triangles.

An example: the traditional linear elements form an 'pyramid' like element around a node. For quadratic elements a degree of freedom corresponding to the value of function on an edge is of second kind. For every degree of freedom the is one linear equation.

Note that the boundary nodes should be treated here as well, although, they are releveant only in the first two cases.

Integration

The integration of the basis elements should be done for only one type of triangle. Use canonical coordinates (see handouts??) and transform the results afterwards. (see codes below).

Helpful Matlab functions

• onboundary.m Tests whether a node, edge or a triangle is on boundary.
• inwhichtri.m Finds a triangle, which has the point inside it. Returns also the canonical coordinates.
• pdeplotn.m Plots the geometry given by pdetool with numbering. (Anybody know if this is already implemented?)

Also following functions in the toolbox directory file:/p/appl/math/matlab/toolbox/pde/ might turn out helpful. E.g.

• pdetrg.m returns the area of each triangle
• pdetridi.m Side lengths and areas of triangles
• pdesubix.m Find indices in index vector
• pdeent.m For a given triangle find indices of the neighbouring triangles.