This repository contains the MATLAB toolbox FEMT2D
.
It provides a discretization of Helmholtz type PDE's in two dimensions using first or second order Lagrange elements on triangular meshes.
Besides a basic MATLAB installation, a mesh generator is required for setting up the triangulation. Currently, only triangular meshes generated by the Triangle library (references see below) are supported.
To make the creation of 2-D triangular meshes more user-friendly, the installation of PyGIMLi is strongly recommended.
FEMT2D
has been tested with MATLAB R2020a. Earlier releases work well also. Also note that compatibility with Octave has not been tested.
To activate rendering of README.md
, download and activate this Chrome extension.
- FEMT2D: FE simulation of 2D Magnetotellurics in MATLAB. https://doi.org/10.5281/zenodo.3955407
- Jonathan Richard Shewchuk, Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator, in "Applied Computational Geometry: Towards Geometric Engineering" (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of Lecture Notes in Computer Science, pages 203-222, Springer-Verlag, Berlin, May 1996.
- Jonathan Richard Shewchuk, Delaunay Refinement Algorithms for Triangular Mesh Generation, Computational Geometry: Theory and Applications 22(1-3):21-74, May 2002.
- Franke, A., Börner, R. U., & Spitzer, K. (2007). Adaptive unstructured grid finite element simulation of two-dimensional magnetotelluric fields for arbitrary surface and seafloor topography. Geophysical Journal International, 171(1), 71-86.
- Börner, R. U. (2010). Numerical modelling in geo-electromagnetics: advances and challenges. Surveys in Geophysics, 31(2), 225-245.
- Rücker, C., Günther, T., Wagner, F.M., 2017. pyGIMLi: An open-source library for modelling and inversion in geophysics, Computers and Geosciences, 109, 106-123, doi: 10.1016/j.cageo.2017.07.011.
Clone this repository.
git clone https://github.com/ruboerner/FEMT2D.git
Then enter the newly created folder FEMT2D
:
cd FEMT2D
For a quick demonstration, simply call
>> demo.driverWeaver
which loads a triangular mesh, and evaluates the response of two quarter-spaces for a given period at predefined observation points. A few typical plots, such as profiles of apparent resistivities, phases, and magnetic transfer functions will be generated.
Inside the FEMT2D
folder, there are a few subfolders. The folders with leading +
, e.g., +fe
, contain packages and have their own namespace.
FEMT2D
├── +demo
├── +fe
├── +mesh
├── +mt
├── +plot
├── +tools
├── Notebooks
├── meshes
From the top-level folder FEMT2D
you can start a demonstration by typing, e.g.,
demo.driverCOMMEMI_2D_0
which will calculate the model responses of model 2D-0 from the COMMEMI model suite.
The folder Notebooks
contains example Jupyter Notebooks which introduce the user into the process of mesh generation using PyGIMLi.
Each of these Notebooks creates a particular set of files in the folder meshes
. The meshes can later be read by routines which reside in the +mesh
folder.
After reading in a Triangle triangulation file, the mesh topology is stored within a MATLAB struct
.
mesh = mesh.getMesh('filename', 'meshes/commemi2d0.1', 'format', 'triangle');
'filename'
corresponds to the basename of a mesh file in the meshes
folder. Note that the actual files have the extensions .poly
, .ele
, and .node
.
If the files in the meshes
folder are, e.g.,
meshes
├── commemi2d0.1.ele
├── commemi2d0.1.node
├── commemi2d0.1.poly
then the corresponding filename
in the parameter list of mesh.getMesh()
would be commemi2d0.1
.
The mesh
structure contains the following fields:
field name | dimension | description |
---|---|---|
node |
2 x np |
node coordinates |
tri2node |
3 x nt |
table of triangles defined by their indices as listed in the array node
|
tri2subdomain |
1 x nt |
array of subdomain number each triangle belongs to |
marker |
1 x nt |
array of attributes associated with nodes |
edge2node |
2 x ne |
table of edges and their node indices |
tri2edge |
3 x nt |
table of edges and their associated triangle index |
nt |
scalar | number of triangles |
np |
scalar | number of nodes |
ne |
scalar | number of edges |
bdEdges |
vector | indices of boundary edges |
bdNodes |
vector | indices of boundary nodes |
Bk |
2 x 2 x nt |
array of matrices defining the affine map |
Binv |
2 x 2 x nt |
array of matrices |
detBk |
nt x 1 |
array of determinants of |
area |
nt x 1 |
array of all area = 0.5 * detBk
|
|
The finite element method associates the triangulation with specific basis functions defined on the elements and identifies the degrees of freedom (DOFS).
freq = 0.01;
sigma = [1e-1, 1e-1, 1e-2, 1e-2, 1e-14];
mu = ones(size(sigma));
fem = fe.FEMproblem('mesh', mesh, ...
'elementtype', 'Lagrange', ...
'order', 2, 'dimension', 2, ...
'sigma', sigma, 'mu', mu, ...
'polarization', 'both', ...
'frequency', freq, ...
'verbose', true);
As one can conclude from the above example, the mesh consists of five subdomains each of which is associated with a unique parameter (here sigma
denotes the electrical conductivity mu
is the relative magnetic permeability
In the given example, the solutions of the Helmholtz equations (for both magnetotelluric polarizations) are desired at a frequency freq
of 0.01 Hz. The polynomial order order
of the Lagrange elements is 2, i.e., we use quadratic (second order) finite element basis functions.
The numerical accuracy is always much better when second-order elements are employed.
After calling FEMproblem
, the fem
structure contains the following fields:
field name | dimension | description | default |
---|---|---|---|
mesh |
struct | copy of mesh struct |
- |
order |
scalar | 1 or 2 for linear or quadratic elements, resp. | 1 |
dimension |
scalar | spatial dimension (currently only 2D is implemented) | 2 |
app |
struct | contains application (PDE) specific parameters, e.g., the frequency frequency for the MT case |
- |
polarization |
string | indicates whether E- ('epol' ), H- ('hpol' ) or both ('both' ) polarizations are desired |
'both' |
elementtype |
type of finite elements | 'Lagrange' |
|
sigma |
nt x 1 |
elementwise electrical conductivities |
- |
mu |
nt x 1 |
elementwise relative permeabilities |
- |
stripair |
scalar | indicator that controls whether the Air halfspace has to be excluded (e.g., for the MT for H-polarization) | - |
istrip |
scalar | subdomain number associated with the Air halfspace that has to be stripped from mesh for MT | - |
coorddofs |
2 x ndofs |
coordinates of the DOFS | - |
elem2dofs |
(order * 3) x nt |
array of elements and their associated DOFS | - |
ndofs |
scalar | number of degrees of freedom (DOFS) | - |
bddofs |
vector | boundary DOF indices | - |
Compl |
ndofs x nbddofs |
mapping from full to boundary DOFS | - |
Null |
ndofs x (ndofs - nbddofs) |
mapping from full to inner DOFS | - |
EnL |
scalar | electric field at |
- |
EnR |
scalar | electric field at |
- |
HnL |
scalar | magnetic field at |
- |
HnR |
scalar | magnetic field at |
- |
The assembly of the FE matrices can be carried out using a call to FEMassemble
:
fem = fe.FEMassemble(fem, 'output', 'matrices', 'verbose', true);
After return, the fem
struct contains the following additional fields:
field name | dimension | description |
---|---|---|
Ke |
ndofs x ndofs |
E-pol stiffness matrix (sparse) |
Me |
ndofs x ndofs |
E-pol mass matrix (sparse) |
Kh |
ndofs x ndofs |
H-pol stiffness matrix (sparse) |
Mh |
ndofs x ndofs |
H-pol mass matrix (sparse) |
Usually the solution (and its gradient) is only required at a few points within the computational domain. To this end we setup observation operators for the desired polarizations and field components, such that
The observation operators can be constructed using the call
yobs = tools.asRow(-99000:1000:99000);
obs = [yobs; 0.01 + zeros(size(xobs))];
nobs = size(obs, 2);
fem = fe.getQ(fem, obs);
The fem
structure now contains an additional field Q
which in turn contains the discrete observation operators.
It is always assumed that the strike direction is aligned with the
field name | dimension | description |
---|---|---|
Qe |
nobs x ndofs |
evaluates |
QeHy |
nobs x ndofs |
evaluates |
QeHz |
nobs x ndofs |
evaluates |
To enable the plot of numerical results across the vertical plane we call getQfull
:
fem = fe.getQfull(fem);
which computes the observation operators necessary for the interpolation of the computed quantities at the triangular elements midpoints. This provides the following additional operators, which are also stored under fem.Q
:
field name | dimension | description |
---|---|---|
Qefull |
nobs x ndofs |
interpolates |
Qhfull |
nobs x ndofs |
interpolates |
QJefull |
nobs x ndofs |
interpolates |
QeHyfull |
nobs x ndofs |
interpolates |
QeHzfull |
nobs x ndofs |
interpolates |
QhEyfull |
nobs x ndofs |
interpolates |
QhEzfull |
nobs x ndofs |
interpolates |
Boundary conditions are required to ensure a unique solution to the PDE. We assume a 1-D layered at both the left and right boundaries. Note that the boundaries may be different!
After successful setup of the FE discretization using FEMproblem
and the assembly of the linear system using FEMassemble
, the Dirichlet boundary values for all four domain boundaries are calculated in the function removeDirichlet
.
A typical call would look like
fem = fe.removeDirichlet(fem, 'sigmaBCL', [1e-2, 1e-1, 1e-3], 'thicknessBCL', [2e4, 2e3], ...
'sigmaBCL', [1e-1, 1e-3], 'thicknessBCL', 2e4, 'frequency', freq);
In the given example, the conductivities and thicknesses at the left and right boundaries are different.
At the top and bottom model boundaries, the Dirichlet values are obtained by horizontal interpolation considering a Neumann-type condition imposed at the four corners of the modelling domain.
If a simple 1-D background with identical layers at the left and the right model boundaries is considered, the two keyword pairs sigmaBCL
and sigmaBCR
as well as thicknessBCL
and thicknessBCR
can be replaced by the keyword pair sigmaBC
and thicknessBC
in the argument list of removeDirichlet
:
fem = fe.removeDirichlet(fem, 'sigmaBC', [1e-2, 1e-3], 'thicknessBC', 2e4, freq);
Although arbitrary topography can be modelled in the computational domain, the Earth's surface at both boundaries is restricted to be at FEMT2D
.
The assembled system of linear equations can be solved using the call
sol = fe.FEMsolve(fem, 'verbose', true);
Depending on the desired polarization, the MATLAB structure sol
contains the degrees of freedom (DOFS) sol.ue
for E-polarization and/or sol.uh
for the H-polarization.
Note that the actual values at the desired points of interest have to be interpolated using the afformentioned Q
operators.
To obtain all possible field components for both polarizations at the positions defined by obs
and one frequency fem.app.frequency
, one would call
iw = -2i * pi * fem.app.frequency
Ex = fem.Q.Qe * sol.ue;
Ey = fem.Q.QhEy * sol.uh;
Ez = fem.Q.QhEz * sol.uh;
Hx = fem.Q.Qh * sol.uh;
Hy = fem.Q.QeHy * sol.ue / iw;
Hz = fem.Q.QeHz * sol.ue / iw;
The last two statements originate from the componentwise evaluation of Maxwell's curl equation
FEMT2D
is able to incorporate inhomogenous magnetic permeabilities, which are taken into account in the forming of fem.Q.QeHy
and fem.Q.QeHz
when calling getQ
.
A plot of the modulus of the electrical current density can be produced using a call to patchplotConst
:
plot.patchplotConst(fem.mesh.node, ...
[fem.elem2dofs(1:3, :); tools.asRow(abs(fem.Q.QJefull * sol.ue))], ...
@(x) x, 'none');
colormap(gca, jet(256));
h = colorbar();
ylabel(h, '|J| in A/m^2');
axis equal tight ij
patchplotConst
takes in its first argument a list of nodes as obtained by triangle. The second argument is an 3 x nt
FE element table augmented with an array of scalar quantities in an additional 4th row.
In our example, this scalar quantitiy is the absolute value of the induced current density,
Examples of function handle | Description |
---|---|
@(x) x |
no transformation is applied, i.e., the data will be displayed as is |
@log10 |
the common logarithm is applied to the data |
The last argument controls the color of the triangular mesh. Colors can be defined either as normalized RGB values, e.g., [0.1, 0.1, 0.1]
or as color abbreviation, e.g., 'w'
for white. The mesh can be suppressed by using 'none'
as mesh color.
The triangular elements are colored according to the color table set by a call to, e.g., colormap(gca, jet(256))
.
All plot routines have been called by the top level function plot.plotMT()
which accepts the following parameters:
Parameter | Value | Purpose |
---|---|---|
fem |
fem struct | |
sol |
sol struct with results of FEMsolve
|
|
mesh |
mesh struct | |
obs |
2 x np array with np observation points |
|
periods |
nf array of periods |
required for sounding curve plots |
profile |
rhoa+phase , tipper
|
plot for a single frequency along profile defined by obs
|
sounding |
rhoa+phase |
sounding curves at a particular station given by its number site
|
section |
Jx |
depth section of a given field component, here |
section |
conductivity |
depth section of conductivities |
section |
subdomain |
depth section of subdomain number which correspond to entries of the sigma and mu arrays required in the call to FEMproblem
|
More plot types will be added to +plot/plotMT.m
in the future.
MIT License
Copyright (C) 2020 Ralph-Uwe Börner
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.