FEMT2D

KO | EN

gitlite — search

#javascript #java #android #php #ruby #python #rails #opengl #shaders #cms #scala #search

★ 45 Open GitHub ↗

FE simulation of 2D Magnetotellurics in MATLAB

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

[![DOI](https://zenodo.org/badge/276864022.svg)](https://zenodo.org/badge/latestdoi/276864022)

# FEMT2D

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.

## Requirements
Besides a basic MATLAB installation, a mesh generator is required for setting up the triangulation. Currently, only triangular meshes generated by the [_Triangle_ library](https://www.cs.cmu.edu/~quake/triangle.html) (references see below) are supported.

To make the creation of 2-D triangular meshes more user-friendly, the installation of [PyGIMLi](http://www.pygimli.org) is strongly recommended.

`FEMT2D` has been tested with MATLAB R2020a. Earlier releases work well also. Also note that compatibility with [Octave](http://www.octave.org) has not been tested.

To activate rendering of $\LaTeX$ expressions within this `README.md`, download and activate this [Chrome extension](https://chrome.google.com/webstore/detail/mathjax-plugin-for-github/ioemnmodlmafdkllaclgeombjnmnbima).

#### References
* 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.

## Installation
Clone this repository.
```bash
git clone https://github.com/ruboerner/FEMT2D.git
```

Then enter the newly created folder `FEMT2D`:
```bash
cd FEMT2D
```

## Quick start

For a quick demonstration, simply call
```matlab
>> 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.

## Folder structure

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.,
```matlab
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.

## Details of the toolbox
After reading in a _Triangle_ triangulation file, the mesh topology is stored within a MATLAB `struct`.

```matlab
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

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 $n_t$ |
| `np`      | scalar | number of nodes $n_p$ |
| `ne`      | scalar | number of edges $n_e$ |
| `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 $B_k$ from reference element to each of the $k=1,\dots,n_t$ domain elements |
| `Binv` | `2 x 2 x nt` | array of matrices $B_k^{-1}$  |
| `detBk` | `nt x 1` | array of determinants of $B_k$, $\text{det}(B_k)$ |
| `area`    | `nt x 1` | array of all $k=1,\dots,n_t$ triangular element surface areas, i.e., `area = 0.5 * detBk` |
|

### The `fem` structure

The finite element method associates the triangulation with specific basis functions defined on the elements and identifies the degrees of freedom (DOFS).

```matlab
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 $\sigma$ given in S/m and `mu` is the relative magnetic permeability $\mu_r$, such that $\mu = \mu_r \mu_0$ with $\mu_0 = 4 \pi \cdot10^{-7}$ Vs/Am).

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 $\sigma$ in S/m | - |
| `mu` | `nt x 1` | elementwise relative permeabilities $\mu_r$ | - |
| `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 $z=0$ for the left boundary | - |
| `EnR`    | scalar | electric field at $z=0$ for the right boundary | - |
| `HnL`    | scalar | magnetic field at $z=0$ for the left boundary | - |
| `HnR`    | scalar | magnetic field at $z=0$ for the right boundary | - |

The assembly of the FE matrices can be carried out using a call to `FEMassemble`:

```matlab
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) |

### The observation operators

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 $\mathbf Q : \mathbb C^N \to \mathbb C^n$ maps from an $N$-dimensional discrete solution (or components of its gradient) to an $n$-dimensional vector of spatial observations. Usually, $N \gg n$ and $\mathbf Q$ is sparse. The observation operators are independent of the frequency.

The observation operators can be constructed using the call
```matlab
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 $x$-axis of a right-handed coordinate system.
More precisely, the $x$-axis points into the plane, the $y$-axis point to the right, and the $z$-axis points downward.

| field name | dimension | description           |
|------------|:----------|-----------------------|
| `Qe`  | `nobs x ndofs` | evaluates $E_x$ for E-pol |
| `QeHy`  | `nobs x ndofs` | evaluates $H_y$ for E-pol |
| `QeHz`  | `nobs x ndofs` | evaluates $H_y$ for E-pol |

To enable the plot of numerical results across the vertical plane we call `getQfull`:
```matlab
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 $E_x$ across the computational domain |
| `Qhfull`  | `nobs x ndofs` | interpolates $H_x$ across the computational domain |
| `QJefull`  | `nobs x ndofs` | interpolates $J_x$ across the computational domain |
| `QeHyfull`  | `nobs x ndofs` | interpolates $H_y$ across the computational domain |
| `QeHzfull`  | `nobs x ndofs` | interpolates $H_z$ across the computational domain |
| `QhEyfull`  | `nobs x ndofs` | interpolates $E_y$ across the computational domain |
| `QhEzfull`  | `nobs x ndofs` | interpolates $E_z$ across the computational domain |

### Boundary conditions

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
```matlab
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`:
```matlab
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 $z=0$ m. This feature will be included in a future release of the `FEMT2D`.

### Solution of the system

The assembled system of linear equations can be solved using the call
```matlab
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 
```matlab
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 
$$\nabla \times \mathbf E = -i\omega\mu \mathbf H.$$
Note that `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`.

### Visualization of the results

A plot of the modulus of the electrical current density can be produced using a call to `patchplotConst`:
```matlab
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, $|J_x|$, for E-polarization.
To transform this quantity, e.g., by using its common logarithm, we apply a function defined by its _function handle_.

| 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 $\|J_x\|$ |
|`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.

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.
```

FEMT2D

Explore Similar Repositories

Was this content helpful?

Recent Feedback

Download README