Finite Element Method with Maple

Copyright 2000, Fedderik van der Bos

Author's note: This worksheet is intended for academic use only. It is not for commerical use.

This worksheet computes solutions of linear second order non-symmetric PDE's using a Finite Element Method (FEM). It also uses the NAG library, which greatly reduces the computing time. Most of the algorithms and also the notation are from the book Introduction to Scientific Computing written by B. Lucquin and O. Pironneau, John Wiley & Sons, 1998. The following problem is used (equation 6.1-3 in the book)

in

on

on

With n the normal on the boundary and is the boundary of .

This worksheet demonstrates use of the accompanying Maple package FEMpackge.mws, which may be downloaded by clicking on the "Download Code" column in the Maple Application Center entry for this application.

The FEM-package contains the following procedures.

• Gridplot - A small procedure for plotting the grid
• Square_Grid - Grid generator on [0,1]x[0,1]
• Generate_Hatfunctions - Hat function generator
• Convert_Data_Interior - Preprocessing inside the domain
• Convert_Data_Interior - Preprocessing on the boundary
• Linear System - Procedure to make the lineair system .

Solving the Linear System can easily be done by LinearSolve from the LinearAlgebra package.

It uses the following global variables

• Vertices - Array for the vertex-coordinates
• Triangles - Array for the triangles
• Boundary - List of several Arrays containing Boundary-segments
• Number_of_Vertices
• Number_of_Triangles
• Number_of_Boundary_segments - List of the number of vertices in the different boundary segments
• M_data, U_data, A0_data, F_data, A1_data, G_data and Phi0_data - Different Arrays in which information is stored during pre-processing. The data is used in building up the linear system.
• Matrix_A and Vector_F - The Matrix and the Vector forming the linear system

The following variables have to be defined.

• From the PDE: M, u, a0, f, a1, g and phi0. a1, g and phi0 have to be defined as an list for every boundary-segment.
• N - Integer controling the number of vertices and triangles.
• Type_Dirichlet - List of booleans (true of false) telling wether a boundary segment has Dirichlet or Robin boundary conditions.
• Bandwidth - Integer for the bandwidth with which the Matrix A is formed. This involves some loss of information, but less memory is needed and calculation can be done faster, in the book see chapter 3. The grid generator Square_Grid(N) is so, that if Bandwidth is set to 2N+1 that no information is lost.
• e - A large real constant for a trick described on page 90 and 91.

Now we can begin.

Libname has to be redefined in order to be able to use the FEM-package on my PC.

> restart;

> libname:="F:\\Program Files\\Maple Libraries",libname:

> with(plots):

> with(linalg):

> with(FEM);

Warning, the name changecoords has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

The usage of Square_Grid, Generate_Hatfunction and Gridplot can be seen in the following section. For the other procedures more data is needed and their usage can been seen in the examples.

Generate Grid and Hatfunctions

Square_Grid and Gridplot

> Square_Grid(5);

> Vertices;

> Triangles;

> Boundary;

> Number_of_Vertices;

> Number_of_Triangles;

> Number_of_Boundary_segments;

> Gridplot();

Generate_Hatfunctions

> Generate_Hatfunctions();

> Hatfunctions[1](x);

> plot3d(Hatfunctions[7](x),x[1]=0..1,x[2]=0..1);

>

Examples

Check from a known function

In this example we make the function solution . Also we define M, u, a0 and a1. f and g are calculated. This is a nice way to see if the program works. As you can see in the example I have chosen the functions are quite random.

Functions

> phi:=x->x[1]^2+sin(2*x[1]*Pi)*x[2]+cos(2*x[2]*Pi)+exp(x[1]*x[2]);

> plot_phi:=plot3d(phi(x),x[1]=0..1,x[2]=0..1,color=blue):

> display(plot_phi,axes=normal);

Interior functions.

> M:=x->matrix(2,2,[1,0,0,1]);

> u:=x->[-x[2],x[1]];

> a0:=x->10*sin(x[1]*Pi/10)*cosh(x[2]*Pi);

Boundary functions, look at the way they are defined.

> phi0:=x->[phi(x),phi(x),phi(x),phi(x)];

> a1:=x->[x[1],x[2],x[1]*x[2],sin(x[2])];

> f:=unapply(-diverge(evalm(M(x)&*grad(phi(x),[x[1],x[2]])),[x[1],x[2]])+innerprod(u(x),grad(phi(x),[x[1],x[2]]))+a0(x)*phi(x),x);

For the calculation of g we need to know the normal.

> N:=[[0,-1],[1,0],[0,1],[-1,0]];

> g:=unapply([seq(evalm(N[k]&*M(x)&*grad(phi(x),[x[1],x[2]]))+a1(x)[k]*phi(x),k=1..4)],x);

Input

> Type_Dirichlet:=[true,false,false,true];

> N:=4;

> Bandwidth:=2*N+1;

> e:=10000;

Grid generation

> Square_Grid(N);

> Generate_Hatfunctions();

Calculation

> Convert_Data_Interior(M,u,a0,f);

> Convert_Data_Boundary(Type_Dirichlet,phi0,a1,g);

Data generated

> M_data;

> U_data;

> A0_data;

> F_data;

> Phi0_data;

> A1_data;

> G_data;

> Linear_System(Type_Dirichlet,Bandwidth,e);

> Matrix_A;

> Vector_F;

> Phi:=convert(LinearAlgebra[LinearSolve](Matrix_A,Vector_F),list);

Plotting

> phi_h:=x->sum(Phi[k]*Hatfunctions[k](x),k=1..Number_of_Vertices);

> plot_phi_h:=plot3d(phi_h(x),x[1]=0..1,x[2]=0..1,color=red,labels=[x,y,z]):

> display({plot_phi_h,plot_phi},axes=normal);

>

Acoustics

This is an nice example from the book (paragraph 1.4). The "+" sign is wrong in the book, it has to be a "-" sign.

Functions

> f:=x->exp(-10*((x[1]-.1)^2+(x[2]-.5)^2));

> g:=x->[0,0,0,0];

> Type_Dirichlet:=[false,false,false,false];

> phi0:=x->0;

> M:=x->matrix(2,2,[1,0,0,1]);

> u:=x->[0,0];

> a0:=x->.8;

> a1:=x->[0,0,0,0];

Input

> N:=5;

> Bandwidth:=2*N+1;

> e:=10000;

Grid generation

We want the grid to look like a concert hall, so the domain [0,1]x[0,1] is transformed.

> Square_Grid(N);

> Trans:=x->[x[1],piecewise(x[1]<=.2,x[2],x[1]<=.8,(.8+x[1])*x[2]+(x[1]-.2),(2.4-x[1])*x[2]+x[1]-.2)];

> VT:=[seq(Trans(Vertices[i]),i=1..Number_of_Vertices)]:

> Vertices:=Array(1..Number_of_Vertices,1..2,VT,datatype=float,readonly);

> Gridplot();

> Generate_Hatfunctions();

Calculation

> Convert_Data_Interior(M,u,a0,f);

> Convert_Data_Boundary(Type_Dirichlet,phi0,a1,g);

> Linear_System(Type_Dirichlet,Bandwidth,e);

> Phi:=convert(LinearAlgebra[LinearSolve](Matrix_A,Vector_F),list);

Plotting

> phi_h:=x->sum(Phi[k]*Hatfunctions[k](x),k=1..Number_of_Vertices);

> plot3d(phi_h(x),x[1]=0..1,x[2]=0..(2.2));

>