The objective of this assignment is to design and implement FEM-based numerical techniques for CEM applications. As part of the assignment, you need to upload to Canvas the report in pdf format, addressing the questions and tasks, described below.
Exercise 1 - FEM in 1D
In this exercise, we will solve numerically the Helmholtz equation we studied in Lecture: FEM 1D (also section in the textbook).
As part of the exercise:
1.1 Describe the advantages and disadvantages of the FEM method with respect to the Finite Difference Method.
1.2 What are the five different steps of deriving the FEM?
1.3 When do we call the FEM method the "Galerkin method"?
1.4 Derive the weighted average of the residual (the equations we solve numerically in the FEM method) by integrating by parts and expanding the solution using nodal basis functions:
From the original Helmholtz problem:
1.5 Implement a code in Matlab, Python or Julia to solve numerically the Helmholtz equation above with
- alpha = 1
- beta = 1
- a = 0
- b = pi
- s(x) = 10*sin(3*x)
You need to perform the following steps (in Matlab code):
1. Define the configuration (number of collocation points, stiffness matrix, ...)
n = 10;
h = pi/n;
% Global stiffness matrix A and right-hand side vector b.
A = zeros(n-1,n-1);
b = zeros(n-1,1);
2. Build A and b iterating over the collocation points:
xim1 = (i-1)*h;
xi = i*h;
% pt1 and pt2 are the points used in 2-point Gauss quadrature for evaluating integrals
pt1 = (xi+xim1)/2 - h/(2*sqrt(3));
pt2 = (xi+xim1)/2 + h/(2*sqrt(3));
% Form Stiffness matrix A
pint = .5*(alpha(pt1)+ alpha(pt2))*h;
qiint = .5*(beta(pt1)*(pt1-xim1)^2 + beta(pt2)*(pt2-xim1)^2)*h;
qim1int = .5*(beta(pt1)*(xi-pt1)^2 + beta(pt2)*(xi-pt2)^2)*h;
if i < n
A(i,i) = A(i,i) + pint/h^2 + qiint/h^2;
if i > 1
A(i-1,i-1) = A(i-1,i-1) + pint/h^2 + qim1int/h^2;
qiim1 = .5*h*(beta(pt1)*(xi-pt1)*(pt1-xim1) + beta(pt2)*(xi-pt2)*(pt2-xim1));
if i < n
A(i-1,i) = -pint/h^2 + qiim1/h^2;
A(i,i-1) = A(i-1,i);
% Right-hand side
fiint = .5*(s(pt1)*(pt1-xim1) + s(pt2)*(pt2-xim1))*h;
fim1int = .5*(s(pt1)*(xi-pt1) + s(pt2)*(xi-pt2))*h;
if i < n
b(i) = b(i) + fiint/h;
if i > 1
b(i-1) = b(i-1) + fim1int/h;
a) In this code our alpha, beta and s depend on x, we can define functions to calculate them:
function alphaval = alpha(x)
alphaval = 1;
function betaval = beta(x)
betaval = 1;
% source term
function sval = s(x)
sin3x = sin(3*x);
sval = 9*sin3x;
sval = sval + sin3x;
b) We solve the integral for preparing the assembly matrix numerically instead of analytically. For this integral, we use a numerical technique, called 2-point Gauss quadrature (an extension of the midpoint trapezoidal rule we saw in the first module).
3. Solve the linear system using Matlab command \ to find the solution f
% Solve the linear system.
f = A\b;
4. Calculate the analytical solution to calculate the error later on
xi = i*h;
ftrue(i,1) = sin(3*xi); % this analytical solution
As part of the exercise:
1.5.1 Report and describe Matlab code to calculate the element of stiffness matrix, Aij, and source term elements bi. Where the derivative of the test/weight function is taken? How do you solve the linear system?
1.5.2 What is the sparsity of the stiffness matrix? Plot the sparsity. Is what you expected?
1.5.3 Perform and plot a convergence test of the method, increasing the number of elements. What is the accuracy of the method?
1.5.4 Change the code to solve the Poisson equation d^2 Phi / dx^2 = - rho in 1D domain [0, pi]. the source term is rho = sin(2x) and Phi(0) = Phi(pi) = 0. Plot the solution of the potential using a different number of elements. It is not necessary to calculate the analytical solution for comparison.
Exercise 2 - FEM2D for Poisson Equation
The goal of this assignment is to use the FEM to solve the same problem we solve in Assignment II (coaxial transmission line problem). To answer the theory questions of this assignment, study Section 6.3 of the textbook.
The assignment is divided into three main parts:
2.1 State and give a brief description of the basic steps for solving the Poisson Equation using the FEM.
2.2 What are the element and the assembly matrices? What is the relation between the element and assembly matrices?
2.3 Download the finite element codes and mesh FEMpoisson.m, CmpElMtx.m and unimesh0.mat.
2.3.1 What is the sparsity of the assembly matrix A? What is its structure?
2.3.2 What are two data structures used for storing information about the mesh? What are the numbers of elements and nodes in the mesh file unimesh0.mat?
2.3.3 Plot the Potential. Hint: check Textbook Section 6.5.6 on how to plot the FEM grid and values.
2.3.4 Measure the error of the FEM code using the capacitance result: you can take the finite difference high-resolution capacitance result and compare it with the FEM capacitance result. What is the Finite-Difference resolution that gives an error that is comparable to the FEM error?
2.3.5 How would the FEM governing equations and solution change if we used symmetry BCs and solved the Poisson equation only in 1/4 of the domain as we did in the Finite-Difference case? You don't need to report the calculations but you can explain how the method would change it and why it would change it.
Bonus Exercise - FEM Method for Electromagnetic Problems - Eigenvalue Problem
In this exercise, we will calculate the eigenmodes H and eigenvalues k2 for a cavity resonator with a circular metal boundary of radius a = 1m. This exercise is discussed in Section 6.5.6 of the textbook.
The solution satisfies the eigenvalue problem nabla x nabla H = k2 H with boundary conditions n x nabla x H = 0.
We will use the following codes: and edgeFEM2D.m (calculating Mass and Stiffness Matrices), plotfield.m (plotting the eigenmodes), mesh_circular_coarse.mat (circular mesh with triangular elements).
3.1 Why are edge elements needed? Why are they called edge elements? Why are they referred to as curl-conforming elements? List some of the characteristic properties of edge elements.
3.2 Which boundary condition corresponds to n x nabla x H = 0?
3.3 What are the numbers of elements and nodes in the mesh mesh_circular_coarse.mat?
3.4 What are Stiffness (S) and Mass (M) matrices when solving an eigenvalue problem with FEM?
3.5 Calculate and plot the first 10 eigenvalues (use edgeFEM2D.m) and compare the results with the results in Fig. 6.28.
3.6 Plot the fundamental modes and the modes associated with the second and third eigenvalues. Comment the results.