Extension Functions

delaunay_points_bnd_rectangle(bmin, bmax)
    -> (Vector{Tuple{Float64,Float64}}, Vector{Vector{Int}})

Create the initial point list and closed outer boundary curve for a rectangular domain, for use as input to triangulate from DelaunayTriangulation.jl.

The arguments bmin = [xmin, ymin] and bmax = [xmax, ymax] give the lower-left and upper-right corners of the bounding rectangle. The four corners are stored as Float64 tuples in counter-clockwise order. The returned bndcurve encodes the closed loop in the format expected by the boundary_nodes keyword of triangulate.

Returns a tuple (points, bndcurve) where:

• points is a Vector{Tuple{Float64,Float64}} containing the four corner vertices,
• bndcurve is a Vector{Vector{Int}} that closes the rectangular boundary.

After this call, constraint edges and interior points can be appended with delaunay_points_add_segment. The result is then passed to triangulate from DelaunayTriangulation.jl, and the triangulation converted with delaunay_to_simplicial to obtain a EuclideanComplex.

This function is only available when DelaunayTriangulation.jl is loaded.

Example

using DelaunayTriangulation
using ConleyDynamics

points, bndcurve = delaunay_points_bnd_rectangle([-5, -5], [5, 5])
tri = triangulate(points; boundary_nodes = bndcurve)
sc  = delaunay_to_simplicial(tri)

delaunay_points_add_nodes(points, newpoints)
    -> Vector{Tuple{Float64,Float64}}

Append a list of interior nodes to the given points for a Delaunay triangulation.

newpoints is a Vector{Vector{<:Real}} where each element [x, y] is a 2D point which describes an interior node. All points in newpoints are appended to points (converted to Float64), and the function returns the updated points.

This function is only available when DelaunayTriangulation.jl is loaded.

Example

using DelaunayTriangulation
using ConleyDynamics

# Rectangular domain with a circular hole and an open cut
points, bndcurve = delaunay_points_bnd_rectangle([-5, -5], [5, 5])

newpoints = [-4 .+ 8 .* rand(2) for k in 1:12]   # 12 random points
points    = delaunay_points_add_nodes(points, newpoints)

tri = triangulate(points; boundary_nodes = bndcurve)
sc  = delaunay_to_simplicial(tri)

delaunay_points_add_segment(points, segments, newseg)
    -> (Vector{Tuple{Float64,Float64}}, Set{Tuple{Int,Int}})
delaunay_points_add_segment(points, newseg)
    -> (Vector{Tuple{Float64,Float64}}, Set{Tuple{Int,Int}})

Append a polygonal curve to the point list and constraint-edge set for a constrained Delaunay triangulation.

newseg is a Vector{Vector{<:Real}} where each element [x, y] is a 2D point on the curve. If the first and last entries of newseg coincide (Euclidean distance less than 1e-9) the curve is treated as closed: only one copy of the shared endpoint is stored and the closing edge connects the last new point back to the first. Otherwise the curve is treated as open and both endpoints are stored separately.

All points in newseg are appended to points (converted to Float64) and all successive pairs are added as directed constraint edges to segments.

The three-argument form takes an existing segments::Set{Tuple{Int,Int}} and extends it. The two-argument form is a convenience wrapper that initialises segments as an empty set; use it when adding the first constraint curve after delaunay_points_bnd_rectangle.

Returns the updated (points, segments).

This function is only available when DelaunayTriangulation.jl is loaded.

Example

using DelaunayTriangulation
using ConleyDynamics

# Rectangular domain with a circular hole and an open cut
points, bndcurve = delaunay_points_bnd_rectangle([-5, -5], [5, 5])

n = 12
theta  = range(0, 2*pi, length = n+1)
circle = [[4.0*cos(t), 4.0*sin(t)] for t in theta]  # closed curve
cut    = [[-2.0, -2.0], [0.0, 1.0], [2.0, -1.0]]    # open curve

points, segments = delaunay_points_add_segment(points, circle)
points, segments = delaunay_points_add_segment(points, segments, cut)

tri = triangulate(points; boundary_nodes = bndcurve, segments = segments)
sc  = delaunay_to_simplicial(tri)

delaunay_points_add_split_segs(points, newsegments; verbose=false)
    -> (Vector{Tuple{Float64,Float64}}, Set{Tuple{Int,Int}})
delaunay_points_add_split_segs(points, segments, newsegments; verbose=false)
    -> (Vector{Tuple{Float64,Float64}}, Set{Tuple{Int,Int}})

Compute a planar straight-line arrangement from newsegments, add all resulting vertices and sub-segments to points and segments, and return the updated pair.

Use this instead of delaunay_points_add_segment when the input curves may intersect each other or have collinear overlaps.

newsegments is a Vector of segments, each given as a length-2 vector of [x, y] endpoints, e.g. [[x1, y1], [x2, y2]]. The function calls split_segments internally to resolve all crossings and overlaps into a non-crossing arrangement, then appends the resulting vertices and edges to the existing points and segments.

The two-argument form is a convenience wrapper that initialises segments as an empty set. The three-argument form takes an existing segments::Set{Tuple{Int,Int}} and extends it.

If the keyword argument verbose is true, the function prints a brief summary: number of input segments, number of intersection points created, and number of output vertices and edges added to points.

Returns the updated (points, segments).

This function is only available when DelaunayTriangulation.jl is loaded.

Example

using DelaunayTriangulation
using ConleyDynamics

points, bndcurve = delaunay_points_bnd_rectangle([-2, -2], [2, 2])

# Two crossing diagonals — split_segments resolves the crossing automatically
segs = [[[-1.0, -1.0], [1.0, 1.0]],
        [[-1.0,  1.0], [1.0, -1.0]]]

points, segments = delaunay_points_add_split_segs(points, segs; verbose=true)

tri = triangulate(points; boundary_nodes = bndcurve, segments = segments)
sc  = delaunay_to_simplicial(tri)

delaunay_to_simplicial(tt; p::Int=2) -> EuclideanComplex

Convert a Delaunay triangulation from DelaunayTriangulation.jl into an EuclideanComplex.

The argument tt is a triangulation object as returned by triangulate from DelaunayTriangulation.jl. Ghost triangles and ghost vertices (the package's boundary sentinel elements) are automatically excluded via the each_solid_* iterators, so the result contains only the geometric simplices.

The keyword argument p specifies the field characteristic of the resulting complex: p=0 gives the rationals, p>0 gives GF(p). The default is p=2.

Vertex labels are zero-padded decimal strings ("001", "002", …) whose width is determined by the total vertex count. The resulting EuclideanComplex carries the embedded 2D coordinates of each simplex.

This function requires DelaunayTriangulation.jl to be loaded before ConleyDynamics.jl.

Example

using DelaunayTriangulation
using ConleyDynamics

pts = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0)]
tri = triangulate(pts)
sc  = delaunay_to_simplicial(tri)
sc3 = delaunay_to_simplicial(tri; p=3)

tetgen_mesh(dom::TetDomain; maxvolume::Union{Real,Nothing}=nothing,
            minratio::Union{Real,Nothing}=1.414,
            mindihedral::Union{Real,Nothing}=nothing,
            steiner_on_boundary::Bool=true, quiet::Bool=true,
            switches::Union{String,Nothing}=nothing) -> TetGenResult

Create a constrained tetrahedral mesh of a domain with TetGen.jl.

The domain dom is described by a TetDomain, see for example tetgen_domain_box and tetgen_domain_add_surface. The keyword arguments control the mesh quality and size:

• maxvolume: Maximum volume of the tetrahedra. If nothing, there is no global volume constraint. Regional volume constraints of the domain are always applied. Note that TetGen treats volume constraints as approximate: typically about one percent of the tetrahedra exceed them, by up to about fifty percent.
• minratio: Maximum ratio of circumradius to shortest edge length of the tetrahedra. Smaller values give better shaped tetrahedra, but more of them. The default 1.414 is the standard value of TetGen. Values below 1 are rejected, since TetGen does not terminate for them. If both minratio and mindihedral are nothing, no quality refinement is performed.
• mindihedral: Minimum dihedral angle of the tetrahedra, in degrees, strictly between 0 and 90. If nothing, only minratio is used. If minratio is nothing, but mindihedral is given, the ratio bound 2.0 is used. Note that TetGen does not guarantee this bound, and that large values can make TetGen fail.
• steiner_on_boundary: If false, no additional points are added on the boundary and constraint facets. Note that this also severely limits the refinement close to these facets.
• quiet: If true, the output of TetGen is suppressed.
• switches: A complete TetGen switch string, which overrides all other keyword arguments. It should contain p if the facets of the domain are to be respected, and it should not contain f, since otherwise interior triangles are reported as boundary triangles with marker 0.

Quality refinement is switched on by default, since flat tetrahedra can make the construction of multivector fields fragile. The values of maxvolume, minratio, and mindihedral have to be finite, otherwise an ArgumentError is raised.

If the facets of the domain intersect, or if the quality bounds are too strict, TetGen fails, and an ArgumentError is raised. In that case, check the constraint surfaces and curves, relax the quality bounds, or use Gmsh.jl for intersecting geometry. The debugging files tetgen-tmpfile* which TetGen writes into the current directory in such cases are removed. An ArgumentError is also raised if the mesh does not contain any tetrahedra, which happens if the outer boundary is not closed, or if a hole point does not lie inside a closed constraint surface.

The resulting TetGenResult can be converted into an EuclideanComplex with tetgen_to_simplicial.

This function requires TetGen.jl to be loaded.

Example

using TetGen, ConleyDynamics

dom = tetgen_domain_box([-2,-2,-2], [2,2,2])
V, T = trisurface_sphere([0,0,0], 1.0, 2)
dom = tetgen_domain_add_surface(dom, V, T)
res = tetgen_mesh(dom; maxvolume = 0.01)
ec  = tetgen_to_simplicial(res)
sphere_cells = tetgen_marker_cells(res, ec, 1)

create_simplicial_delaunay3d(bmin::AbstractVector{<:Real},
                             bmax::AbstractVector{<:Real},
                             pdist::Real, attempts::Int;
                             p::Int=2, box_corners::Bool=false,
                             rng::AbstractRNG=Random.default_rng())
                             -> EuclideanComplex

Create a three-dimensional Delaunay triangulation of a random point sample inside a box. The complex is over the rationals if p=0, and over GF(p) if p>0.

The function selects a random sample of points inside the box [bmin[1],bmax[1]] x [bmin[2],bmax[2]] x [bmin[3],bmax[3]], with minimum distance pdist between the points, see poisson_points_box for the meaning of attempts. It then computes the Delaunay tetrahedralization of these points with TetGen.jl, and returns it as an EuclideanComplex. The random number generator can be specified via rng for reproducible results.

As in the planar case create_simplicial_delaunay, the complex covers only the convex hull of the sample, so it does not fill the box completely near the boundary, and the tetrahedra on the boundary of the convex hull can be very flat. If box_corners=true, the eight corners of the box are added to the sample, so that the convex hull is the whole box. For quality meshes of a box, use tetgen_mesh with tetgen_domain_box instead.

This function requires TetGen.jl to be loaded.

TetDomain

Description of a three-dimensional domain for tetrahedral mesh generation with TetGen.jl, as a piecewise linear complex.

The struct has the following fields:

• points::Vector{NTuple{3,Float64}}: The vertices of the domain description.
• facets::Vector{Vector{Int}}: The facets, each given as a vector of indices into points. A facet with at least three vertices is a planar polygon, a facet with two vertices is a segment, and a facet with one vertex is a single point. All of them are contained in the final mesh.
• markers::Vector{Int}: One marker for every facet. By convention, the outer boundary has marker 0, and the k-th added constraint surface or curve has marker k.
• holes::Vector{NTuple{3,Float64}}: For every hole, one point inside it. The region around this point which is enclosed by facets is not meshed.
• regions::Vector{Tuple{NTuple{3,Float64},Int,Float64}}: For every region, a point inside it, an integer attribute, and the maximum tetrahedron volume in that region. A nonpositive maximum volume means that there is no regional volume constraint.

A TetDomain is usually created with tetgen_domain_box or tetgen_domain_from_surface, and then extended with the functions tetgen_domain_add_nodes, tetgen_domain_add_surface, tetgen_domain_add_curve, tetgen_domain_add_hole, and tetgen_domain_add_region. None of these functions require TetGen.jl. The domain is then meshed with tetgen_mesh.

All facets have to be non-intersecting, except along shared vertices and edges. In particular, constraint surfaces cannot cross each other or the outer boundary.

TetGenResult

Tetrahedral mesh created by tetgen_mesh.

The struct has the following fields:

• points::Vector{Vector{Float64}}: The vertex coordinates of the mesh, including all points added by the mesh generator.
• tetrahedra::Vector{Vector{Int}}: The tetrahedra, as vectors of four indices into points.
• trifaces::Vector{Vector{Int}}: The boundary and constraint triangles, as vectors of three indices into points.
• trifacemarkers::Vector{Int}: The marker of the facet which contains each triangle in trifaces.
• tetattributes::Vector{Int}: The region attribute of every tetrahedron, or an empty vector if the domain has no regions.
• curvesegments::Vector{Tuple{NTuple{3,Float64},NTuple{3,Float64},Int}}: The segments of all constraint curves of the domain, given by their endpoints and their marker.

The mesh can be converted into an EuclideanComplex with tetgen_to_simplicial, and the cells on constraint surfaces and curves, or in regions, can be found with tetgen_marker_cells and tetgen_region_cells.

tetgen_domain_box(bmin::AbstractVector{<:Real}, bmax::AbstractVector{<:Real};
                  markers::AbstractVector{<:Integer}=zeros(Int,6)) -> TetDomain

Create a TetDomain for the box [bmin[1],bmax[1]] x [bmin[2],bmax[2]] x [bmin[3],bmax[3]].

The domain consists of the eight corners of the box and its six faces, which are quadrilateral facets with outward orientation. By default, all faces have marker 0. The optional argument markers assigns individual markers to the faces, in the order x = bmin[1], x = bmax[1], y = bmin[2], y = bmax[2], z = bmin[3], and z = bmax[3]. This makes it possible to find the cells on some of the faces with tetgen_marker_cells, for example to describe an exit set. This is the three-dimensional analogue of delaunay_points_bnd_rectangle.

tetgen_domain_from_surface(V::AbstractVector, T::AbstractVector) -> TetDomain

Create a TetDomain whose outer boundary is the closed triangulated surface given by the vertices V and the triangles T.

The surface is usually created by one of the trisurface_* functions, and it has to be closed, i.e., every edge has to belong to exactly two triangles. Otherwise an ArgumentError is raised. Vertices with the same coordinates are identified before this check, so the surface can also be given as a list of separate triangles, each with its own copies of the vertices, as in STL files. Vertices which do not belong to any triangle are ignored. For example, the output of trisurface_sphere gives a ball, and the output of trisurface_torus gives a solid torus. All facets get marker 0.

tetgen_domain_add_nodes(dom::TetDomain, newpoints::AbstractVector) -> TetDomain

Add nodes to a TetDomain.

The vector newpoints contains the coordinates of the new nodes, as vectors or tuples of length three. Each of these nodes inside the domain becomes a vertex of the final mesh. Nodes outside of the domain or in holes are ignored by the mesh generator, and nodes which coincide with points of the domain are merged with them. Nodes must not lie on a facet or a constraint curve without being one of its vertices, since TetGen reports such nodes as self-intersections and fails. The function returns a new TetDomain, and the original domain is not changed.

tetgen_domain_add_surface(dom::TetDomain, V::AbstractVector, T::AbstractVector;
                          marker::Union{Int,Nothing}=nothing,
                          hole::Bool=false,
                          region::Union{Tuple{Int,Real},Nothing}=nothing)
                          -> TetDomain

Add a triangulated constraint surface to a TetDomain.

The surface is given by its vertices V and its triangles T, for example as created by one of the trisurface_* functions. Its triangles become facets of the domain, so the final mesh will contain the surface as a union of triangles. The facets get the marker marker, and if no marker is given, the next free marker maximum(dom.markers) + 1 is used, or 1 if this maximum is zero or less. The cells of the final complex on the surface can be found with tetgen_marker_cells.

For closed surfaces, two additional options are available, which cannot be combined:

• hole=true removes the inside of the surface from the domain.
• region=(attr, maxvol) marks the inside of the surface as a region with integer attribute attr, and with maximum tetrahedron volume maxvol (use a nonpositive value for no volume constraint). The tetrahedra in this region can be found with tetgen_region_cells.

In both cases, the surface has to be closed and consistently oriented, and a point inside the surface is determined automatically. This point avoids the surfaces which are already part of the domain, so that a region bounded by a sphere around a previously added hole is the shell between the two surfaces. If the surface consists of several connected components, such as two disjoint spheres, one point is determined for every component, so that all of them become holes or regions. Note that for nested components, such as two concentric spheres, both the inner ball and the shell between the spheres are removed or marked. If this is not intended, add the surfaces separately, starting with the inner one, or use tetgen_domain_add_hole or tetgen_domain_add_region with explicit points.

Vertices of the surface which coincide exactly with each other or with existing points of the domain are identified, before the surface is checked. Triangles which then have repeated vertices are rejected, and vertices which do not belong to any triangle are ignored. Apart from such shared vertices and edges, the surface must not intersect other surfaces or the boundary of the domain. For intersecting surfaces, use Gmsh.jl instead, see the manual. The function returns a new TetDomain, and the original domain is not changed.

tetgen_domain_add_curve(dom::TetDomain, curve::AbstractVector;
                        marker::Union{Int,Nothing}=nothing) -> TetDomain

Add a polygonal constraint curve to a TetDomain.

The vector curve contains the consecutive points of the curve, as vectors or tuples of length three. If the first and the last point agree up to 1e-9 times the size of the curve, the curve is closed. Every segment of the curve becomes a facet of the domain, so the final mesh will contain the curve as a union of edges. Repeated consecutive points and segments which occur more than once are ignored. The facets get the marker marker, and if no marker is given, the next free marker is used as in tetgen_domain_add_surface.

Vertices of the curve which coincide exactly with existing points of the domain are identified with them. Apart from such shared vertices, the curve must not intersect any facets of the domain, and in particular it cannot lie inside a face of the outer boundary. Curves outside of the domain or inside holes are ignored by TetGen, and tetgen_mesh prints a warning in this case. The function returns a new TetDomain, and the original domain is not changed.

tetgen_domain_add_hole(dom::TetDomain, point::AbstractVector{<:Real}) -> TetDomain

Add a hole to a TetDomain.

The point point has to lie inside a region of the domain which is completely enclosed by facets. This region is then removed from the mesh. A hole point which is not enclosed by facets removes the whole domain. The function returns a new TetDomain.

tetgen_domain_add_region(dom::TetDomain, point::AbstractVector{<:Real},
                         attr::Int, maxvol::Real=-1.0) -> TetDomain

Add a region to a TetDomain.

The region is the part of the domain which contains point and which is bounded by facets. All tetrahedra in this region get the integer attribute attr, and they can be found with tetgen_region_cells. If maxvol is positive, the tetrahedra in the region will have volume at most maxvol, up to the tolerance of TetGen described in tetgen_mesh. The function returns a new TetDomain.

tetgen_to_simplicial(res::TetGenResult; p::Int=2) -> EuclideanComplex

Convert a tetrahedral mesh created by tetgen_mesh into an EuclideanComplex. The complex is over the rationals if p=0, and over GF(p) if p>0.

Vertices of the mesh which do not belong to any tetrahedron are removed. Vertex labels are zero-padded decimal strings as in tetmesh_to_euclidean.

tetgen_marker_cells(res::TetGenResult, ec::EuclideanComplex, marker::Int;
                    closure::Bool=true) -> Vector{Int}

Find the cells of a complex which lie on the facets with a given marker.

The complex ec has to be created from res via tetgen_to_simplicial, otherwise an ArgumentError is raised. The function returns the indices of all triangles of ec which lie on facets of the domain with marker marker, as well as all edges on constraint curves with this marker. If closure=true (default), all faces of these cells are included as well, so the result is a closed subcomplex. For a constraint surface, this gives the surface as a two-dimensional subcomplex, and for the marker 0 it gives the outer boundary of the domain.

The returned cells can be used, for example, to define exit sets or to restrict a multivector field to a constraint surface.

tetgen_region_cells(res::TetGenResult, ec::EuclideanComplex, attr::Int;
                    closure::Bool=true) -> Vector{Int}

Find the tetrahedra of a complex which lie in the regions with a given attribute.

The complex ec has to be created from res via tetgen_to_simplicial, otherwise an ArgumentError is raised, and the domain has to contain regions, see tetgen_domain_add_region and tetgen_domain_add_surface. The function returns the indices of all tetrahedra with region attribute attr. If closure=true (default), all faces of these tetrahedra are included as well, so the result is a closed subcomplex.

Note that TetGen assigns new attributes, larger than all given ones, to the parts of the domain which do not contain a region point. For example, if a single region with attribute 1 is defined inside a closed constraint surface, then the tetrahedra outside of this surface get attribute 2. The attributes of all tetrahedra are stored in res.tetattributes.

gmsh_to_simplicial(; dim::Int=3, p::Int=2,
                   physical::Union{Tuple{Integer,Integer},Nothing}=nothing)
                   -> EuclideanComplex

Convert the mesh of the current Gmsh.jl model into an EuclideanComplex. The complex is over the rationals if p=0, and over GF(p) if p>0.

The function reads the mesh elements of dimension dim from the current model, which has to be meshed already, for example by calling gmsh.model.mesh.generate(3). For dim=3, these are the tetrahedra, and the result is a tetrahedral mesh. For dim=2, the triangles are used, which gives a surface complex embedded in three dimensions, and for dim=1, the edges are used, which gives a graph. All faces of these elements are added automatically. Since Gmsh only stores the triangles on the surfaces and the edges on the curves of the model, and not the interior faces of the tetrahedra, dim=2 gives only the triangles on the surfaces of the model, and dim=1 only the edges on its curves. If physical=(pdim, ptag) is given, only the elements of the physical group with dimension pdim and tag ptag are used, and pdim has to equal dim.

Only linear simplices are supported, i.e., edges, triangles, and tetrahedra whose nodes are their vertices. If the elements of dimension dim include other types, an ArgumentError is raised. This happens for quadrangles, hexahedra, prisms, and pyramids, which Gmsh creates with its recombination options. It also happens for second order elements, which gmsh.model.mesh.setOrder(2) creates by adding nodes on the edges to follow curved geometry. Only the elements of dimension dim are checked, so that, for example, a hexahedral mesh can still be converted with dim=1. Nodes which do not belong to any of the used elements are dropped. Vertex labels are zero-padded decimal strings as in tetmesh_to_euclidean.

If two distinct nodes of the used elements have the same coordinates, a warning is printed. This usually means that touching or intersecting entities were meshed separately, without gmsh_common_refinement or a gmsh.model.occ.fragment operation, so that the complex is not connected across them. A warning is also printed if some of the used entities of dimension dim, i.e., surfaces for dim=2 and volumes for dim=3, contain no mesh elements. This happens if the meshing failed for them, so that the complex has holes there.

The function only reads the current model. It neither initializes nor finalizes Gmsh, so the model can be used afterwards, for example to find the cells of physical groups with gmsh_marker_cells.

This function requires Gmsh.jl to be loaded.

Example

using Gmsh: gmsh
using ConleyDynamics

gmsh.initialize()
gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.2)
gmsh.model.mesh.generate(3)
ec = gmsh_to_simplicial()
gmsh.finalize()

gmsh_file_to_simplicial(filename::AbstractString; dim::Int=3, p::Int=2)
                        -> EuclideanComplex

Read a mesh or geometry file with Gmsh.jl and convert it into an EuclideanComplex. The complex is over the rationals if p=0, and over GF(p) if p>0.

The file can be in any format which Gmsh can open, such as a mesh file (.msh) or a geometry file (.geo). If the file does not contain mesh elements of dimension dim, a mesh of this dimension is generated with the mesh options specified in the file. The mesh is then converted with gmsh_to_simplicial.

If Gmsh is not initialized, the function initializes it and finalizes it at the end. Otherwise, the file is opened in the running Gmsh session. This creates a new model, which becomes the current model, while the previous models remain available.

This function requires Gmsh.jl to be loaded.

gmsh_marker_cells(ec::EuclideanComplex, dim::Integer, tag::Integer;
                  physical::Bool=true, closure::Bool=true) -> Vector{Int}

Find the cells of a complex which belong to a physical group or an entity of the current Gmsh.jl model.

The complex ec has to be created from the current model, for example with gmsh_to_simplicial. If physical=true (default), the function returns the indices of the cells of ec formed by the mesh elements of the physical group with dimension dim and tag tag. If physical=false, the elements of the model entity with dimension dim and tag tag are used instead. Depending on dim, these are vertices, edges, triangles, or tetrahedra. Elements which are not cells of ec are ignored. If closure=true (default), all faces of these cells are included as well, so the result is a closed subcomplex.

Physical groups are the most convenient way to keep track of geometry after boolean operations such as gmsh.model.occ.fragment, since these operations renumber the entities. The cells of ec are identified by the coordinates of their vertices, so the numbering of the nodes does not matter.

This function requires Gmsh.jl to be loaded.

gmsh_common_refinement(objects::AbstractVector;
                       names::Union{AbstractVector{<:AbstractString},Nothing}=nothing)
                       -> Vector{Int}

Split the given geometric objects of the current Gmsh.jl model along all their mutual intersections, and create one physical group for every object.

Every element of objects describes one geometric object by a single entity (dim, tag) or by a vector of such entities, all of the same dimension. For example, a box is given by (3, box), a rectangle by (2, rect), and the surface of a sphere by the vector of its boundary surfaces, see the manual. The objects have to be created with the OpenCASCADE kernel gmsh.model.occ.

The function calls gmsh.model.occ.fragment for all entities of all objects, which replaces them by pieces that only meet along common faces, edges, and points. Afterwards, the mesh of every object consists of the mesh elements of its pieces, and any two objects share exactly the mesh elements of their intersection. In other words, the mesh is a common refinement of all objects. For every object, the function then creates a physical group of the same dimension, which contains all pieces of the object, and it returns the tags of these groups in the order of objects. If names is given, it contains a name for every physical group.

After meshing, the cells of the objects are found with gmsh_marker_cells, and the cells on the intersection of two objects are the intersection of their cell lists. Pieces of an object which lie outside of all volumes, for example the part of a plane outside of a box, are not part of the tetrahedral mesh, and their elements are ignored by gmsh_marker_cells.

This function requires Gmsh.jl to be loaded.

Example

using Gmsh: gmsh
using ConleyDynamics

gmsh.initialize()
box    = gmsh.model.occ.addBox(-2, -2, -2, 4, 4, 4)
sphere = gmsh.model.occ.addSphere(0, 0, 0, 1.0)
plane  = gmsh.model.occ.addRectangle(-2, -2, 0.3, 4, 4)
gmsh.model.occ.synchronize()
sphsurf = gmsh.model.getBoundary([(3, sphere)], false, false)
gmsh.model.occ.remove([(3, sphere)])

pgbox, pgplane, pgsphere = gmsh_common_refinement([(3, box), (2, plane), sphsurf])
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.25)
gmsh.model.mesh.generate(3)

ec     = gmsh_to_simplicial()
circle = intersect(gmsh_marker_cells(ec, 2, pgplane), gmsh_marker_cells(ec, 2, pgsphere))
gmsh.finalize()

gmsh_add_trisurface(V::AbstractVector, T::AbstractVector)
                    -> Vector{Tuple{Int,Int}}

Add a triangulated surface to the current Gmsh.jl model, so that it can be used as an object of gmsh_common_refinement.

The surface is given by the vertices V and the triangles T, for example as created by one of the trisurface_* functions, or as read from an STL file. Every triangle becomes a flat surface of the OpenCASCADE kernel gmsh.model.occ, and adjacent triangles share their points and edges, so the surface is the piecewise linear interpolation of its vertices. The function returns the entities (2, tag) of these surfaces, which form one object of gmsh_common_refinement. The surface does not have to be closed or consistently oriented, and an edge may belong to more than two triangles.

Vertices with the same coordinates are identified, so the surface can also be given as a list of separate triangles, each with its own copies of the vertices. Vertices which do not belong to any triangle are ignored. An ArgumentError is raised for triangles with repeated vertices, for degenerate triangles whose vertices lie on a line, for triangles which occur more than once, and for triangles with an edge of length at most 2e-7, since OpenCASCADE identifies points closer than 1e-7. If OpenCASCADE fails nevertheless, the entities created so far are removed from the model before the error is passed on.

Since every edge of the triangulation becomes a curve of the model, the mesh near the surface is at least as fine as the triangulation, even if Mesh.MeshSizeMax is larger. The triangulation should therefore not be finer than needed. The function calls gmsh.model.occ.synchronize().

This function requires Gmsh.jl to be loaded.

Example

using Gmsh: gmsh
using ConleyDynamics

gmsh.initialize()
box  = gmsh.model.occ.addBox(-2, -2, -2, 4, 4, 4)
V, T = trisurface_parametric((x, y) -> [x, y, (x^2 + y^2) / 2 - 1],
                             (-1.5, 1.5), (-1.5, 1.5), 12, 12)
surf = gmsh_add_trisurface(V, T)

pgbox, pgsurf = gmsh_common_refinement([(3, box), surf])
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.25)
gmsh.model.mesh.generate(3)

ec    = gmsh_to_simplicial()
cells = gmsh_marker_cells(ec, 2, pgsurf)
gmsh.finalize()