FEATool
v1.8 Finite Element Analysis Toolbox |

Grid

After a model geometry has been defined, a computational grid or mesh must be generated to allow for the finite element discretization. This section describes how to create or import a suitable grid.

Grid mode can be selected by pressing the mode button or corresponding menu option. In grid mode the toolbar buttons allow for grid generation, refinement, and accessing the advanced grid generation settings.

- calls the selected grid generation function to generate and unstructured. The desired overall mean
**Grid Size**is specified in the corresponding edit field which can be set manually or with the slider control. - The button is used to uniformly refine a grid. This will split triangular and quadrilateral (2D) grid cells into four subcells, and tetrahedral and hexahedral (3D) ones into eight subcells each. Grid points on boundaries will also be aligned with any geometries in 2D.
- Advanced grid settings can be controlled by pressing the button which opens the
**Grid Generation Settings**dialog box

The grid settings allows selecting the **Grid Generation Algorithm**, either the built-in *DistMesh* algorithm, or the external grid generators *Triangle* (only in 2D) and *Gmsh*. The **Subdomain Grid Size** and **Boundary Grid Size** edit fields can be used to prescribe maximum grid sizes for subdomains and boundaries. The **Quality** slider allows setting the target mesh quality (supported by *DistMesh* and *Triangle*). The number of grid **Smoothing** steps can be specified in the corresponding edit field. **Cell Type** allows selecting between unstructured triangular or structured quadrilateral grid cells (supported by *DistMesh* and *Gmsh*). And lastly, the **Optional Parameters** fields allows for an additional string of property/value pairs to be supplied to the grid generation call.

The **Grid** menu options allows for accessing and specifying the **Grid Generation Call...** which is a MATLAB/Octave function call to support external grid generation tools instead of the default one. **Convert Grid Cells** is used to convert between triangular and quadrilateral cells in 2D, and between tetrahedral and hexahedral cells in 3D.

The **Grid** menu also supports the following import and export of grids through external ASCII format files

**Import Grid > FEniCS/Dolfin Format...**imports a grid in FEniCS/Dolfin format (.xml file) (calls impexp_dolfin).**Import Grid > GiD Format...**imports a grid in GiD format (.msh file) (calls impexp_gid).**Import Grid > Gmsh Format...**imports a grid in Gmsh format (.msh file) (calls impexp_gmsh).**Import Grid > GMV Format...**imports a grid in General Mesh Viewer (GMV) format (.gmv file) (calls impexp_gmv).**Import Grid > Triangle Format...**imports a 2D Triangle grid (.node/.ele files) (calls impexp_triangle).

**Export Grid > FEniCS/Dolfin Format...**exports a grid in FEniCS/Dolfin format (.xml file) (calls impexp_dolfin).**Export Grid > GiD Format...**exports a grid in GiD format (.msh file) (calls impexp_gid).**Export Grid > Gmsh Format...**exports a grid in Gmsh format (.msh file) (calls impexp_gmsh).**Export Grid > GMV Format...**exports a grid in General Mesh Viewer (GMV) format (.gmv file) (calls impexp_gmv).**Export Grid > OpenFOAM Format...**exports a grid in OpenFOAM format (calls grid2foam).**Export Grid > Triangle Format...**exports a 2D Triangle polygonal geometry description (.poly file) (calls impexp_triangle).

Note that in one dimension (1D) the only available option is to create a line grid.

FEATool Multiphysics comes with built-in support for the external mesh generators Gmsh and Triangle. Using either *Gmsh* or *Triangle* typically allows for faster and more robust grid generation, and also supports specifying different grid sizes on subdomains and boundaries. The *Grid Generation Algorithm* can be selected with the corresponding option in the **Grid Generation Settings** dialog box.

If FEATool cannot find the Gmsh or Triangle binaries, FEATool will try to automatically download, (compile), and install them. To install them manually, download *Gmsh* from http://gmsh.info/, and *Triangle* from http://www.cs.cmu.edu/~quake/triangle.html, then extract (compile) and install the binaries into the

lib/gmsh/ lib/triangle/

subdirectories of the FEATool installation directory.

- Gmsh and Triangle can be used in the FEATool GUI as well as on the MATLAB and Octave command lines with the gridgen_gmsh and gridgen_triangle commands.
- External grid generators works by exporting the current FEATool geometry in its own native format, make a system call to the grid generator binary, and then re-import the generated grid.
- For geometries with multiple and overlapping geometry objects the actual subdomain numbering will generally not correspond to the geometry object numbering (two intersecting geometry objects will for example create three or more subdomains and several internal boundaries). In this case the actual subdomain and boundary numbering for vector valued
*hmax*and*hmaxb*arrays can easiest be visualized and determined by first creating a coarse grid and switching to*Equation/Subdomain*and*Boundary*modes, respectively. *Gmsh*propagates the*hmax*and*hmaxb*values down to the specific nodes in the mesh which means that it is currently not possible to exactly define mesh sizes for subdomains and boundaries.

This section describes the format of the grid data structure that FEATool employs as well as advanced command line (CLI) usage such as manually creating and importing grids.

The grid format used by FEATool is specified in the **grid** struct with the following fields

Field | Description | Size |
---|---|---|

p | Grid point coordinates | (n_sdim, n_points) |

c | Grid cell connectivity | (n_edges_per_cell, n_cells) |

a | Grid cell adjacency | (n_edges_per_cell, n_cells) |

b | Boundary information | (3+n_sdim, n_boundary_points) |

s | Grid cell subdomain list | (1, n_cells) |

The coordinates of the grid points are specified by an array **p**, with the row number indicating the *i*-th coordinate direction, and column number the corresponding grid point number.

Cell connectivities are specified in the **c** array, where each column point to the grid points (in **p**) making up each cell. The row index gives the local vertex number and the column index the cell number. Moreover, the grid points are numbered counterclockwise in each cell.

Adjacency, meaning pointers to neighboring cells, are specified in the **a** array. Similar to **c** the row index gives the local edge or face number (starting with the corresponding grid point in **c**) and the column index points to the cell number. If the edge is on an external boundary the corresponding entry in **a** is *0*.

Boundary information is specified in the **b** array. The cell, edge/face, and boundary numbers are prescribed in the first to third rows. The last *n_sdim* rows consists of the outward pointing unit normals corresponding to each boundary cell edge (or face in 3D).

A list of subdomain numbers for each cell is given in the **s** vector.

The various grid cells available in FEATool are defined in this sub section.

In one dimension only the simple straight line segment grid cell is available. Grid Example 1 shows how this can be defined and used.

The 2D triangular grid cell is defined with vertices in counter clockwise order. Local edges e_{i} are defined by the local vertices v_{i} as

[ e1 ; [ v1 v2 ; e2 ; = v2 v3 ; e3 ] v3 v1 ];

The quadrilateral grid cell is similarly defined with vertices in counter clockwise order. Local edges e_{i} are defined by the local vertices v_{i} as

[ e1 ; [ v1 v2 ; e2 ; v2 v3 ; e3 ] = v3 v4 ; e4 ] v4 v1 ];

The 3D tetrahedral grid cell is defined with vertices in counter clockwise order. Local edges e_{i} are defined by the local vertices v_{i} as

[ e1 ; [ v1 v2 ; e2 ; v2 v3 ; e3 ; v3 v1 ; e4 ; = v1 v4 ; e5 ; v2 v4 ; e6 ] v3 v4 ];

The local faces f_{i} are defined by the local vertices v_{i} as

[ f1 ; [ v1 v2 v3 ; f2 ; v1 v2 v4 ; f3 ; = v2 v3 v4 ; f4 ] v3 v1 v4 ];

The hexahedral brick grid cell is defined with vertices in counter clockwise order. Local edges e_{i} are defined by the local vertices v_{i} as

[ e1 ; [ v1 v2 ; e2 ; v2 v3 ; e3 ; v3 v4 ; e4 ; v4 v1 ; e5 ; v1 v5 ; e6 ; v2 v6 ; e7 ; = v3 v7 ; e8 ; v4 v8 ; e9 ; v5 v6 ; e10 ; v6 v7 ; e11 ; v7 v8 ; e12 ] v8 v1 ];

The local faces f_{i} are defined by the local vertices v_{i} as

[ f1 ; [ v1 v2 v3 v4 ; f2 ; v1 v2 v6 v5 ; f3 ; v2 v3 v7 v6 ; f4 ; = v3 v4 v8 v7 ; f5 ; v4 v1 v5 v8 ; f6 ] v5 v6 v7 v8 ];

The following code can be used to manually define a one dimensional grid with 10 uniformly spaced cells on the line (0..1)

grid.p = 0:0.1:1; grid.c = [1:10;2:11]; grid.a = [0:9;2:10 0]; grid.b = [1 10;1 2;1 2;-1 1]; grid.s = ones(1,10);

Alternatively one can use the grid utility function linegrid

grid = linegrid( 10, 0, 1 );

The plotgrid function can be used to plot and visualize a grid

plotgrid( grid )

A 2 by 2 unit square rectangular grid can be created with the following commands

grid.p = [repmat([0,0.5,1],1,3);0 0 0 0.5 0.5 0.5 1 1 1]; grid.c = [1 2 5 4;2 3 6 5;4 5 8 7;5 6 9 8]'; grid.a = [0 2 3 0;0 0 4 1;1 4 0 0;2 0 0 3]; grid.b = [1 2 2 4 4 3 3 1;1 1 2 2 3 3 4 4;1 1 2 2 3 3 4 4; ... 0 0 1 1 0 0 -1 -1;-1 -1 0 0 1 1 0 0]; grid.s = ones(1,4);

The rectgrid function can also be used to create rectangular grids, in this case

grid = rectgrid( 2, 2, [0 1;0 1] );

As before the grid can be plotted, visualizing both grid point and cell numbers, with

plotgrid( grid, 'nodelabels', 'on', 'cellabels', 'on' )

Similarly, the boundaries can be visualized with the function plotbdr (subdomains can be visualized with plotsubd)

plotbdr( grid )

As FEATool also supports simplex triangular cells in two dimensions a grid consisting of quadrilaterals can easily be converted with the utility function quad2tri

grid = quad2tri( grid )

Reverse conversions can be made with tri2quad function which uses Catmull-Clark subdivision. In 3D the hex2tet and tet2hex functions perform similar conversions.

A more complex example, a grid for a rectangle with a circular hole can be created by first creating geometry objects (a rectangle and a circle), applying a formula to construct the geometry shape, and then calling the automatic unstructured grid generation function gridgen

% Geometry definition. xmin = 0; xmax = 1; ymin = 0; ymax = 1; tag1 = 'R1'; gobj1 = gobj_rectangle( xmin, xmax, ymin, ymax, tag1 ); xc = (xmax-xmin)/2; yc = (ymax-ymin)/2; r = 0.25; tag2 = 'C1'; gobj2 = gobj_circle( [xc yc], r, tag2 ); geom.objects = { gobj1 gobj2 }; formula = 'R1 - C1'; geom = geom_apply_formula( geom, formula ); fea.geom = geom; % Grid generation. hmax = 0.1; fea.grid = gridgen( fea, 'hmax', hmax );

As before the grid can be plotted, visualizing both grid point and cell numbers, with

plotgrid( grid, 'nodelabels', 'on', 'cellabels', 'on' )

FEATool allows importing and exporting finite element grids between FEniCS, GiD, Gmsh, General Mesh Viewer (GMV), OpenFOAM, and Triangle formats with the corresponding impexp_dolfin, impexp_gid, impexp_gmsh, impexp_gmv, impexp_foam, and impexp_triangle functions. These functions can also be accessed from the **Grid** and **Postprocessing** menus of the FEATool GUI.

Due to the simple grid format syntax it is possible to manually import grids from other software. The process essentially consists of first exporting the grid point coordinates and grid cell connectivity data from the external grid generation tool into separate text files. Then import them into Octave or MATLAB, after which they can be reshaped and used by FEATool. The following steps describe the process

- The first step is to output the grid to an ASCII text format. If possible save the grid output to two files, one for the grid vertex coordinates, and another with the grid cell connectivities (this specifies which grid points/vertices make up each cell). If this is not possible one will have to manually open the grid output file in a text editor, and cut and save the grid coordinates and cell connectivities to two different files.
Load the two files in Octave/MATLAB (here it is assumed they are saved as

*p.txt*and*c.txt*):load p.txt load c.txt

- Reshape the grid coordinates (
**p**variable) so that it has the form where is the number of space dimensions (1, 2, or 3) and is the total number of grid points (**p**essentially consists of rows of x, y, and z-coordinates). Precisely how to reshape depends on the output format from the external grid generator, one might not have to do anything (check the shape by entering the command**size(p)**or**whos**), it might be enough to transpose as**p = p';**, or one might have to reshape with something like**p = reshape(p,n_sdim,n_p);**. - Similarly, reshape the cell connectivity array
**c**so that it has the shape where is the number of vertices on each cell (for example 3 for triangles) and is the total number of cells. Each column should point to the corresponding grid points in**p**that make up the cell. - The grid vertices must run in counterclockwise order on each cell, so reorient them if necessary. This is usually accomplished by changing the row order, for example
**c = c([3 2 1],:);**which the ordering for triangles. Use the function gridadj to create an array that points to neighboring cells for each cell edge

n_sdim = size(p,1); a = gridadj( c, n_sdim );

Create a vector that assigns a subdomain number for each cell. If the geometry should be one single block and thus not split, a unit row vector is sufficient

n_c = size(c,2); s = ones( 1, n_c );

Use the function gridbdr to create boundary edge/face information. The

*gridbdr*function splits external boundaries where the angles are acuteb = gridbdr( p, c, a );

Alternatively, one can manually construct and set the boundary numbering, for example as

[ix_edge_face,ix_cells] = find( a==0 ); b = [ix_cells'; ix_edge_face'; ones(1,length(ix_cells))]; b = [b; gridnormals(p,c,b)];

this creates a boundary array

**b**where all external boundaries are joined. With additional boundary information one can now split the boundary edges and faces (columns) by assigning different boundary numbers (third row in**b**). A helpful function to use here is findbdr which allows one to find boundary numbers and column indices in to**b**by prescribing logical expressions such as*x>0.5*. For example[~,ix] = findbdr( struct('p',p,'c',c,'b',b), 'x>0.5', false ); b(3,ix) = 2;

sets boundary number

*2*on all boundary edges/faces where*x>0.5*is true.Create a grid

*struct*to hold all the grid informationgrid.p = p; grid.c = c; grid.a = a; grid.s = s; grid.b = b;

- Now the grid can be used by FEATool functions and subroutines on the command line.
Optionally, one can also import the grid into the FEATool GUI by using the Menu option

File > Import > Variables From Main Workspace...

Select the created**grid**variable to import and press**Import**. This loads it into the local FEATool memory workspace and can now be accessed from the FEATool Terminal (the bottom command **>>** edit field in the main GUI window). By entering the commandfea.grid = grid;

in the FEATool Terminal command line (not the main Octave/MATLAB command window) effectively replaces the current grid with the imported one. Press the**Grid**mode button to update and show the new grid.

The computational finite element library in FEATool supports FEM shape functions for structured grids (quadrilaterals in 2D and hexahedra in 3D). Although more difficult to generate automatically, structured grids are often computationally superior, allowing for higher accuracy with a smaller number of cells.

In FEATool, structured tensor-product grids of basic geometric shapes can easily be generated on the command line with the following Octave/MATLAB m-script functions

Function | Description |
---|---|

linegrid | Create a 1D line grid |

circgrid | Create a 2D structured circular grid of quadrilateral cells |

holegrid | Create a 2D rectangular grid with a circular hole |

rectgrid | Create a 2D rectangular grid of quadrilateral cells |

ringgrid | Create a 2D grid of a ring shaped domain |

blockgrid | Create a 3D structured block grid |

cylgrid | Create a 3D structured cylindrical grid |

spheregrid | Create a 3D grid for a spherical domain |

Moreover, the grid utility functions selcells, delcells, gridextrude, gridmerge, gridrevolve, gridrotate, and gridscale can be used to manually modify, transform, and join grids to more complex structures.

The figure below shows three examples of more complex grids created by using these functions.

The first flow over cylinder benchmark grid is for example created with the following commands:

grid1 = ringgrid( [0.05 0.06 0.08 0.11 0.15], ... 32, [], [], [0.2;0.2] ); grid2 = holegrid( 8, 1, [0 0.41;0 0.41], 0.15, [0.2;0.2] ); grid2 = gridmerge( grid1, 5:8, grid2, 1:4 ); grid3 = rectgrid( [0.41 0.5 0.7 1 1.4 1.8 2.2], ... 8, [0.41 2.2;0 0.41] ); grid = gridmerge( grid3, 4, grid2, 6 );

And the lower right revolved grid with these commands:

grid = rectgrid(); grid = gridscale( grid, {'1-(y>0.5).*(y-0.5)' 1} ); grid = delcells( selcells( ... '((x<=0.8).*(x>=0.2)).*(y<=0.2)', grid ), grid ); grid = gridrevolve( grid, 20, 0, 1/4, 2, pi/2, 0 );

The last example with two brackets attached to an I-beam section is more complex:

grid01 = ringgrid( 1, 20, 0.03, 0.06, [0;0] ); indc01 = selcells( grid01, 'y<=sqrt(eps)' ); grid01 = delcells( grid01, indc01 ); grid02 = holegrid( 5, 1, .06*[-1 1;-1 1], .03, [0;0] ); indc02 = selcells( grid02, 'y>=-sqrt(eps)' ); grid02 = delcells( grid02, indc02 ); grid2d = gridmerge( grid01, [5 6], grid02, [7 8] ); grid1 = gridextrude( grid2d, 1, 0.02 ); grid1 = gridrotate( grid1, pi/2, 1 ); grid2 = grid1; grid1.p(2,:) = grid1.p(2,:) + 0.03; grid2.p(2,:) = grid2.p(2,:) - 0.01; x_coord = [ -0.08 linspace(-0.06,0.06,6) 0.08 ]; y_coord = [ -0.2 -0.15 -0.1 -0.05 -0.03 -0.01 ... 0.01 0.03 0.05 0.1 0.15 0.2 ]; grid3 = blockgrid( x_coord, y_coord, 1, ... [-0.08 0.08;-0.2 0.2;-0.08 -0.06] ); grid4 = blockgrid( 1, y_coord, 5, ... [-0.01 0.01;-0.2 0.2;-0.18 -0.08] ); grid5 = grid3; grid5.p(3,:) = grid5.p(3,:) - 0.12; grid = gridmerge( grid1, 9, grid3, 6 ); grid = gridmerge( grid2, 9, grid, 17 ); grid = gridmerge( grid4, 6, grid, 23, 1 ); grid = gridmerge( grid5, 6, grid, 29, 2 );

After, the grids have been created on the command line they can also be imported into the FEATool GUI (by using the **File > Import > Variables From Main Workspace option**, after which the command **fea.grid = grid;** needs to be entered at the FEATool command prompt).

Using quadrilateral grid cells are often advantageous to simplex and triangular cells in that they can provide somewhat more accuracy when aligned with geometry features and also tends to require less overall grid cells. Quadrilateral grid generation can be specified by selecting *Quad* grid **Cell Type** option in the Grid Settings dialog box (available with *DistMesh* and *Gmsh* in 2D). Quadrilaterals unfortunately do not tend to allow for easy automatic grid generation although It is also possible to subdivide triangles and convert into quads the resulting grids are often of poor quality.

The optional gridgen_quad command line function was originally designed to align quadrilateral cell edges with immersed interfaces described by distance and level set functions (analog to DistMesh). The algorithm then uses the zero level set contour from the distance functions to align grid cell edges with external geometry object boundaries. Furthermore, the cells are split in a way to treat edge cases such as when an interface segment crosses a cell diagonal.

The following functions are available for working with and processing grids

Function | Description |
---|---|

gridgen | Unstructured automatic grid simplex generation (DistMesh) |

gridgen_gmsh | External grid generation with Gmsh (2D and 3D) |

gridgen_triangle | External grid generation with Triangle (2D) |

gridgen_quad | Unstructured automatic grid quadrilateral generation |

gridrefine | Refine a grid uniformly |

gridsmooth | Apply smoothing to a grid |

gridextrude | Extrude a grid from 2D to 3D |

gridrevolve | Extrude and revolve a 2D grid to 3D |

gridrotate | Rotate a grid along a specified axis |

gridscale | Apply scaling to a grid |

gridmerge | Merge two grids along boundaries |

quad2tri | Convert a grid of quadrilateral cells to triangular cells |

tri2quad | Convert a grid of triangular cells to quadrilateral cells |

hex2tet | Convert a grid of hexahedral cells to tetrahedral cells |

tet2hex | Convert a grid of tetrahedral cells to hexahedral cells |

impexp_dolfin | Import and export FEniCS/Dolfin grid data |

impexp_foam | Import and export OpenFOAM grid data |

impexp_gid | Import and export GiD grid data |

impexp_gmv | Import and export GMV grid and postprocessing data |

impexp_gmsh | Import Gmsh grid and postprocessing data |

impexp_triangle | Import and export 2D Triangle grid and polygon data |

gridcheck | Check grid for errors |

gridadj | Create grid adjacency information |

gridbdr | Create grid boundary information |

gridbdrx | Reconstruct internal boundaries |

gridvert | Create grid vertex information |

gridedge | Create grid edge information |

gridface | Create grid face information |

findbdr | Find boundary indices and numbers |

assign_bdr | Sort boundaries according to geometry objects |

project_bdrp | Project boundary points to geometric boundary |

selcells | Find cell indices from an expression |

delcells | Delete a selection of cells from a grid |

plotbdr | Plot and visualize boundaries |

plotsubd | Plot and visualize subdomains |

plotgrid | Plot and visualize grid |