FE simulation of 2D Magnetotellurics in MATLAB.zip

[](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 DO