The Meshutils module

Module containing functionality to read and process tetrahedral meshes in gmsh or nastran format.

Definition of the mesh type

Fields

• name::String: the name of the mesh.
• points::Array: 3×N array containing the coordinates of the N points defining the mesh.
• lines::List: List of simplices defining the edges of the mesh
• triangles::List: List of simplices defining the surface triangles of the mesh
• int_triangles::List: List of simplices defining the interior triangles of the mesh
• tetrahedra::List: List of simplices defining the tetrahedra of the mesh
• domains::Dict: Dictionary defining the domains of the mesh. See comments below.
• file::String: path to the file containing the mesh.
• tri2tet::Array: Array of length length(tetrahedra) containing the indices of the connected tetrahedra.
• dos: special field meant to contain symmetry information of highly symmetric meshes.

Notes

The meshes are supposed to be tetrahedral. All simplices (lines, triangles, and tetrahedra) are stored as lists of simplices. Simplices are lists of integers containing the indices of the points (i.e. the column number in the points array) forming the simplex. This means a line is a two-entry list, a triangle a three-entry list, and a tetrahedron a four-entry list. For convenience certain entities of the mesh can be further defined in the domains dictionary. Each key defines a domain and maps to another dictionary. This second-level dictionary contains at least two keys: "dimension" mapping to the dimension of the specified domain (1,2, or 3) and "simplices" containing a list of integers mapping into the respective simplex lists. More keys may be added to the dictionary to define additional and/or custom information on the domain. For instance the compute_size! function adds an entry with the domain size.

mesh=Mesh(file_name::String; scale=1, inttris=false)

read a tetrahedral mesh from gmsh or nastran file into mesh. The optional scaling factor scale may be used to scale the units of the mesh. The optional parameter inttris toggles whether triangles are sortet into two seperate lists for surface and internal triangles.

collect_lines!(mesh::Mesh)

Populate the list of lines in mesh.lines.

data,surf_keys,vol_keys=color_domains(mesh::Mesh,domains=[])

Create data fields containing an integer number corresponding to the local domain.

Arguments

• mesh::Mesh: mesh to be colored
• domains::List: (optional) list of domains that are to be colored. If empty all domains will be colored.

Returns

• data::Dict: Dictionary containing a key for each colored domain plus magic keys "__all_surfaces__" and "__all_volumes__" holding all colors in one field. These magic fields may not work correctly if the domains are not disjoint.
• surf_keys::Dict: Dictionary mapping the surface domain names to their indeces.
• vol_keys::Dict: Dictionary mapping the volume domain names to their indeces.

Notes

The datavariable is designed to flawlessly work with vtk_write for visualization with paraview. Check the "Interpret Values as Categories" box in paraview's color map editor for optimal visualization.

See also: vtk_write

compute_size!(mesh::Mesh,dom::String)

Compute the size of the domain dom and store it in the field "size" of the domain definition. The size will be a length, area, or volume depending on the dimension of the domain.

full_mesh=extend_mesh(mesh::Mesh, doms; sym_name="Symmetry", blch_name="Bloch", unit=false)

Create full mesh or unit cell from half cell represented in mesh.

Arguments

• mesh::Mesh: mesh representing the half cell. The mesh must span a sector of 2π/2N, where N is the (integer) degree of symmetry of the full mesh.
• doms::List: list of 2-tuples containing the domain names and their copy_degree (see notes below).
• sym_name::String="Symmetry": name of the domain of mesh that forms the symmetry plane of the half-mesh.
• blch_name::String="Bloch": name of the domain of mesh that forms the remaining azimuthal plane of the half-mesh.
• unit::Bool=false: toggle whether extend the mesh to unit cell only.

Returns

• full_mesh::Mesh: representation of the full mesh.

Notes

The routine copies only domains that are specified in doms. These domains are extended according to the specified copy_degree. The following are available:

• :full: extent the domain and save it under the same name in full_mesh.
• :unit: extent the domain and save the individual unit cells labeled from 0 to N-1 in full_mesh.
• :half: extent the domain and save the individual unit cells as half cells labeled from 0 to N-1 in full_mesh where one half-cell contains _img in its name.

Multiple styles can be mixed in one domain specification. An example for doms would be doms=[("Interior", :full), ("Outlet", :unit), ("Flame", :half)]

tet_idx=find_tetrahedron_containing_point(mesh::Mesh,point)

Find the tetrahedron containing the point point and return the index of the tetrahedron as tet_idx. If the point lies on the interface of two or more tetrahedra, the returned tet_idx will be the lowest index of these tetrahedra, i.e. the index depends on the ordering of the tetrahedra in the mesh.tetrahedra. If no tetrahedron in the mesh encloses the point tet_idx==0.

generate_field(mesh::Mesh,func;el_type=0)

Generate field from function funcfor mesh mesh. The element type is either el_type=0 for field values associated with the mesh tetrahedra or el_type=1 for field values associated with the mesh vertices. The function must accept three input arguments corresponding to the three space dimensions.

normal_vectors=get_normal_vectors(mesh::Mesh,output::Bool=true)

Compute a 3×length(mesh.triangles) Array containing the outward pointing normal vectors of each of the surface triangles of the mesh mesh. The vectors are not normalised but their norm is twice the area of the corresponding triangle. The optional parameter output toggles whether a progressbar should be shown or not.

surface_points, tri_mask, tet_mask = get_surface_points(mesh::Mesh,output=true)

Get a list surface_points of all point indices that are on the surface of the mesh mesh. The corresponding lists tri_mask andtet_maskcontain lists of indices of all triangles and tetrahedra that are connected to the surface point. The optional parameteroutput` toggles whether a progressbar should be shown or not.

link_triangles_to_tetrahedra!(mesh::Mesh)

Find the tetrahedra that are connected to the triangles in mesh.triangles (and mesh.inttriangles) and store this information in mesh.tri2tet.

new_mesh=octosplit(mesh::Mesh)

Subdivide each tetrahedron in the mesh mesh into 8 tetrahedra by splitting each edge at its center.

Notes

The algorithm introduces length(mesh.lines) new vertices. This yields a finer mesh featuring 8*size(mesh,2) tetrahedra, 4*length(mesh.triangles) triangles, and 2*length(mesh.lines) lines. From the 3 possible subdivision of a tetrahedron the algorithm automatically chooses the one that minimizes the edge lengths of the new tetrahedra. The point labeling of new_mesh is consistent with the point labeling in mesh, i.e., the first size(mesh.points,2) points in new_mesh.points are identical to the points in mesh.points. Hence, mesh and new_mesh form a hierachy of meshes.

unify!(mesh::Mesh ,new_domain::String,domains...)

Create union of domains and name it new_domain. domains must share the same dimension and new_domain must be a new name otherwise errors will occur.

vtk_write(file_name, mesh, data)

Write vtk-file containing the datavalues data given on the usntructured grid mesh.

Arguments

• file_name::String: name given to the written files. The name will be preceeded by a specifier. See Notes below.
• mesh::Mesh: mesh associated with the data.
• data::Dict: Dictionary containing the data.

Notes

The routine automatically sorts the data according to its type and writes it in up to three diffrent files. Data that is constant on a tetrahedron goes into "$(filename)_const.vtu", data that is linearly interpolated on a tetrahedron goes into "$(filename)_lin.vtu", and data that is quadratically interpolated on a tetrahedron goes into "$(filename)_quad.vtu"

Simple struct element to store additional information on symmetry for rotational symmetric meshes.

#Fields

• DOS::Int64: degree of symmetry
• naxis::Int64: number of grid points lying on the center axis
• nxbloch::Int64: number of grid points lying on the Bloch plane but not on the center axis
• nbody::Int64: number of interior points in half cell
• shiftbody::Int64: difference from body point to reflected body point
• nxsymmetry::Int64: number of grid points on symmetry plane (without center axis)
• nxsector::Int64: number of grid points belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
• naxis_ln::Int64: number of line segments lying on the center axis
• nxbloch_ln::Int64: number of line segments belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
• nxsector_ln::Int64: number of line segments belonging to a unit cell (without cenetraxis, Bloch and Bloch image plane)
• nxsector_tri::Int64: number of surface triangles belonging to a unit cell (without Bloch and Bloch image plane)
• nxsector_tet::Int64: number of tetrahedra of a unit cell
• n: unit axis vector of the center axis
• `pnt: foot point of the center axis
• unit::Bool: true if mesh represents only a unit cell

cmpr=compare(A,B)

Compare two simplices following lexicographic rules. If B has lower rank, then A cmpr ==-1. If B has higher rank than A, cmpr==1. If A and B are identical, cmpr==0.

R = create_rotation_matrix_around_axis(n,α)

Compute the rotation matrix around the directional vector nrotating by an angle α (in rad).

foot = find_foot_of_perpendicular(pnt,pln)

Compute the position of the foot of the perpendicular from the point pnt to the plane pln.

Notes

Code adapted from https://www.geeksforgeeks.org/mirror-of-a-point-through-a-3-d-plane/

p,n = find_intersection_of_two_planes(pln1,pln2)

Find a parametrization of the axis defined by the intersection of the planes pln1and pln2. The axis is parameterized by a point p and direction n.

Notes

The algorithm follows John Krumm's solution for finding a common point on two intersecting planes as it is explained in https://math.stackexchange.com/questions/475953/how-to-calculate-the-intersection-of-two-planes To simplify the code, here the point p0 is the origin.

idx = find_smplx(list,smplx)

Find index of simplex in ordered list of simplices. If the simplex is not contained in the list, the returned index is nothing

ln_idx = get_line_idx(ln,naxis_ln,nxsector_ln)

Get index ln_idx of line ln in a mesh that features naxis grid points and naxis_ln' lines on the center line,nxsectorgrid points andnxsector_ln' lines in a unit cell excluding the center axis.

ln_idx = get_line_idx(mesh::Mesh, ln)

Get index ln_idx of line ln in mesh mesh.

s = get_line_sector(ln)

Get index s of sector containing line ln, in a mesh that features naxis grid points on the center line and nxsector grid points in a unit cell excluding the center axis..

sidx = get_line_sector(mesh::Mesh, ln)

Get sector idx sidxof line ln in mesh mesh.

sector = get_point_sector(pnt_idx,naxis, nxsector)

Get the sector containing point with index pnt_idx, in a mesh that features naxis grid points on the center line and nxsector grid points in a unit cell excluding the center axis.

sidx=get_point_sector(mesh::Mesh,pnt_idx)

Get sector idx sidxof point with index pnt_idxin mesh mesh.

ridx = get_reflected_index(idx,naxis,nxbloch,nbody,shiftbody,nxsymmetry)

Get index ridx of a point in the first unit cell created from reflection of the point with index idx. The definig mesh feature naxis points on the center axis, nxbloch points on the bloch plane (excluding the center axis), nbody points in the interiror of the first half-cell, nxsymmetrypoints on the symmetry plane (excluding the center axis) nxsector points in a unit cell excluding the center axis. The difference of the indeces of a refernce body point and a reflected body point is shiftbody. The number of points per unit cell (excluding the center axis) is nxsector.

ridx = get_rotated_index(idx,sector, naxis, nxsector, DOS)

Get index ridx of point in sector sector created from rotation of the point with index idx. The definig mesh feature naxis points on the center axis, nxsector points in a unit cell excluding the center axis, and has a degree of symmetry DOS.

insert_smplx!(list,smplx)

Insert simplex in ordered list of simplices.

check = is_point_in_plane(pnt, pln)

Check whether point "pnt" is in plane . Note that there is currently no handling of round off errors.

makenormaloutwards!(pln,testpoint)

Ensure that the parameterization of the plane "pln" has a normal that is directed twoards testpoint. If necessary reparametrize pln.

read a mesh from ansys cfx file

points, lines, triangles, tetrahedra, domains=read_msh2(file_name)

Read gmsh's .msh-format version 2. This is an old format and wherever possible version 4 should be used.

http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format

points, lines, triangles, tetrahedra, domains=read_msh4(file_name)

Read gmsh's .msh-format version 4.

http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format

points, lines, triangles, tetrahedra, domains=read_nastran(filename)

Read simplicial nastran mesh.

Note

There is currently no sanity check for non-simplicial elements. These elements are just skipped while reading.

p = reflect_point_at_plane(pnt,pln)

Reflect point pnt at plane pln and return as p.

Notes

Code adapted from https://www.geeksforgeeks.org/mirror-of-a-point-through-a-3-d-plane/

sort_smplx(list,smplx)

Helper function for sorting simplices in lexicographic manner. Utilize divide an conquer strategy to find a simplex in ordered list of simpleces.

pln=three_points_to_plane(A)

Compute equation for a plane from the three points defined in A. The plane is parameterized as a*x+b*y+c*z+d==0 with the coefficients returned as pln=[a,b,c,d].

Arguments

• A::3×3-Array : Array containing the three points defining the plane as columns.

Notes

Code adapted from https://www.geeksforgeeks.org/program-to-find-equation-of-a-plane-passing-through-3-points/