Creating Regular Grids with Triangulate.jl and Extendable Grids

I am trying to solve a PDE system on an unstructured mesh. I currently use the Triangulate.jl to create the mesh and then use OrdinaryDiffEq.jl to solve the system of ODEs, one for each cell of the mesh. Now I am interested in seeing if I get any advantage by solving the system on an unstructured mesh(of triangles) over using a structured grid (of rectangles). I have already looked in to ExtendableGrids.jl which says it can do it. However, its documentation is a bit difficult to comprehend. Does anyone know how to generate a structured grid or a grid using only equilateral triangles, using these packages?

Now I am interested in seeing if I get any advantage by solving the system on an unstructured mesh(of triangles) over using a structured grid (of rectangles).

ExtendableGrids.jl in the moment only supports meshes of triangles (even in the case of rectangular geometries). We have plans to extend this to more general cases. So both via Triangulate and ExtendableGrids you will get meshes of triangles.

Could you elaborate a bit more on your goals ?

BTW, for a PDE solver working on these grids together with DifferentialEquations.jl, you may have a look VoronoiFVMDiffEq.jl.

PS: Yeah, I am aware, the documentation is not optimal in the moment.

Thank you for your reply! I am currently looking at your lectures to use VoronoiFVM.jl. I will also have a look at VoronoiFVMDiffEq.jl afterwards. However, they don’t directly address my goals, which are:

  1. Solve the PDE on a 2d rectangular mesh.
  2. Solve the PDE on a 2d mesh of equilateral triangles (which can also be substituted by the example here, as I just need the orthogonality of edges to the line joining the cell centres). Also, from your comment on GitHub, I think this can be done using the angle(-q) switch of Triangulate.jl and SimplexGridFactory.jl.

As I have stated before, I have already solved the system using unstructured meshes from Trianglulate.jl now I want to see the advantage that I gain by using an unstructured mesh in a domain with non-rectangluar boundary. For that I need to solve the problem on a structured grid (grids mentioned in points 1, 2) using the same tools and compare the results (accuracy, computation time).

PS: An example of how to use the angle switch in Trianglulate.jl would be very helpful.

I just need the orthogonality of edges to the line joining the cell centres

Ah this would just come from the the Delaunay property.
Just to be precise, please distinguish:

  • mesh of equilateral triangles (all angles 60°)
  • mesh of acute triangles (all angles <90°)
  • mesh of weakly acute triangles (all angles <=90°)
  • Delaunay mesh: (sum of angles opposite to a given edge <=180°)
  • Boundary conforming Delaunay mesh: Delaunay + all circumcenters in the closure of the computational domain.

Equilateral, acute and weakly acute are all boundary conforming Delaunay.

There is no widely acknowledged software for creating acute/weakly acute triangulations. OTOH algorithms for Delaunay are well understood (in 2D), and Triangle (behind Triangulate.jl) is one implementation.

My questions:

  • Which of these types of meshes do you really need ?
  • Where do you place the unknowns of your discretization ? triangle nodes or triangle (circum?)centers ?

(EDITs: fix angle error, and Delaunay allows for <=…)

I need the First one all angles 60 degrees and Boundary Conforming Delaunay. As I am placing the unknown on the circum centres, so I require them to be inside the domain.

Please realize that not all domains can be covered with equilateral triangles. What you are asking is impossible in general domains.

Yes, I now realize the impossibility of the grid that I was imagining. Thank you!

Thanks for understanding that I meant 60° :sweat_smile:

1 Like

No issues! But according to the documentation of triangle the ‘q’ switch may not converge for angles larger than 33°. So can provide the interim grid points via Triangulate.pointlist be a suitable alternative?

Yes, this should work. It is not necessary that all points should be part of some segment.

I recently implemented a triangular structured mesher for my PSSFSS code. The triangular mesh is structured in that it is derived from a plaid rectangular mesh by adding a diagonal to each rectangle. The code is pretty simple and would be easy to adapt to your needs, I believe. Here is the documentation string for the function:

"""
    make_plaid_mesh(xr, yr, area, ntri, is_inside) -> sheet::RWGSheet
Generate a structured, plaid triangular mesh from list of required coordinates and predicate function
# Input Arguments
- `xr`, `yr`: Vectors of required x and y coordinates for vertices of the geometry to be meshed.
- `area`:  The area of the geometry to be meshed.
- `ntri`:  The desired number of triangles for the area to be meshed.
- `is_inside`: A predicate function where `is_inside(x,y) -> tf::Bool` determines whether a point (x,y)
  is within the region to be meshed.
#  Return value:
- `sheet`: A variable of type RWGSheet with fields ρ, e1, e2, fe, and fv properly initialized. The
  mesh results from a plaid rectangular tesselation containing at least the vertices in the Cartesian
  product of `xr` and `yr`, the latter supplemented with additional points to refine the mesh, and then
  converted to a triangular tesselation by adding a diagonal to each rectangle.
"""

If it looks useful to you, you can find the code here. Note that the file structuredtri.jl containing this function is included in the file Elements.jl which contains the necessary using statements.

1 Like