Application Center - Maplesoft

App Preview:

Finite element methods for solving PDEs

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

asymmetricPDE.mws

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)

-diverge(M*grad(phi))+innerprod(u,grad(phi))+a[0]*p... in Omega

Phi = Phi[0] on Gamma[1]

innerprod(M*grad(phi),n)+a[1]*phi = g on Gamma[2]

With n the normal on the boundary and Gamma = Gamma[1]+Gamma[2] is the boundary of Omega .

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 A*Phi = F .

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

[Convert_Data_Boundary, Convert_Data_Interior, Gene...

>

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

`Number of Vertices`, 61, `Number of Triangles`, 10...

> Vertices;

_rtable[1647548]

> Triangles;

_rtable[31567820]

> Boundary;

_rtable[31526492]

> Number_of_Vertices;

61

> Number_of_Triangles;

100

> Number_of_Boundary_segments;

[6, 6, 6, 6]

> Gridplot();

[Maple Plot]

Generate_Hatfunctions

> Generate_Hatfunctions();

61, `Hatfunctions generated`

> Hatfunctions[1](x);

PIECEWISE([1.000000000-5.000000000*x[1]-5.000000000...
PIECEWISE([1.000000000-5.000000000*x[1]-5.000000000...
PIECEWISE([1.000000000-5.000000000*x[1]-5.000000000...
PIECEWISE([1.000000000-5.000000000*x[1]-5.000000000...
PIECEWISE([1.000000000-5.000000000*x[1]-5.000000000...

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

[Maple Plot]

>

Examples

Check from a known function

In this example we make the function solution phi . 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]);

phi := proc (x) options operator, arrow; x[1]^2+sin...

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

> display(plot_phi,axes=normal);

[Maple Plot]

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

M := proc (x) options operator, arrow; matrix(2,2,[...

u := proc (x) options operator, arrow; [-x[2], x[1]...

a0 := proc (x) options operator, arrow; 10*sin(1/10...

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])];

phi0 := proc (x) options operator, arrow; [phi(x), ...

a1 := proc (x) options operator, arrow; [x[1], 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);

f := proc (x) options operator, arrow; -2+4*sin(2*x...
f := proc (x) options operator, arrow; -2+4*sin(2*x...

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

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

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

g := proc (x) options operator, arrow; [-sin(2*x[1]...
g := proc (x) options operator, arrow; [-sin(2*x[1]...
g := proc (x) options operator, arrow; [-sin(2*x[1]...
g := proc (x) options operator, arrow; [-sin(2*x[1]...

Input

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

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

> N:=4;

> Bandwidth:=2*N+1;

> e:=10000;

N := 4

Bandwith := 9

e := 10000

Grid generation

> Square_Grid(N);

`Number of Vertices`, 41, `Number of Triangles`, 64...

> Generate_Hatfunctions();

41, `Hatfunctions generated`

Calculation

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

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

`Interior Data ready`

`Boundary Data ready`

Data generated

> M_data;

_rtable[31585700]

> U_data;

_rtable[31588500]

> A0_data;

_rtable[31527412]

> F_data;

_rtable[31527532]

> Phi0_data;

_rtable[31527612]

> A1_data;

_rtable[31528332]
_rtable[31528332]

> G_data;

_rtable[31529092]
_rtable[31529092]

> Linear_System(Type_Dirichlet,Bandwidth,e);

`Linear System ready`

> Matrix_A;

_rtable[31637988]

> Vector_F;

_rtable[31529852]

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

Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...
Phi := [2.00017365524517832, 2.06286050187951408, 1...

Plotting

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

phi_h := proc (x) options operator, arrow; sum(Phi[...

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

[Maple Plot]

>

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];

f := proc (x) options operator, arrow; exp(-10*(x[1...

g := proc (x) options operator, arrow; [0, 0, 0, 0]...

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

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];

phi0 := 0

M := proc (x) options operator, arrow; matrix(2,2,[...

u := proc (x) options operator, arrow; [0, 0] end p...

a0 := .8

a1 := proc (x) options operator, arrow; [0, 0, 0, 0...

Input

> N:=5;

> Bandwidth:=2*N+1;

> e:=10000;

N := 5

Bandwith := 11

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

`Number of Vertices`, 61, `Number of Triangles`, 10...

> 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)];

Trans := proc (x) options operator, arrow; [x[1], p...

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

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

Vertices := _rtable[1491004]

> Gridplot();

[Maple Plot]

> Generate_Hatfunctions();

61, `Hatfunctions generated`

Calculation

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

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

> Linear_System(Type_Dirichlet,Bandwidth,e);

`Interior Data ready`

`Boundary Data ready`

`Linear System ready`

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

Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....
Phi := [.266080080852142020, .265616676732855772, ....

Plotting

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

phi_h := proc (x) options operator, arrow; sum(Phi[...

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

[Maple Plot]

>