Hi there,
I’m trying to automate the discritization of FE meshes for High voltage cables embedded in the ground. This is to conduct thermal analysis and parameterise the geometry of the cable (cables with tri cores).
I wanted to move to Julia, but I think the meshing facilities are not mature enough unfortunately. I’ve tried to use this package called meshes.jl which seems promising. However, I’m finding it difficult to even create the meshing within the cable (an offset circle within a bigger circle. Does anyone know what I’m doing wrong?
using Meshes
import CairoMakie as Mke
Y = [3 * cos(π/20 * (j-1)) for j in 1:40];
X = [3 * sin(π/20 * (j-1)) for j in 1:40];
A = Array{Tuple{Float64, Float64}}(undef, 0)
for i in 1:40
Tup = (X[i],Y[i])
push!(A,(Tup))
end
Y = [ 1+cos(π/20 * (j-1)) for j in 1:40];
X = [ sin(π/20 * (j-1)) for j in 1:40];
B = Array{Tuple{Float64, Float64}}(undef, 0)
for i in 1:40
Tup = (X[i],Y[i])
push!(B,(Tup))
end
ring = Ring(A)
ring1 = Ring(B)
Poly = PolyArea([A,B])
boundary = [A, B]
discretize(Poly ) |> viz
sorry Juliohm , i am fairly new to this and lack knowledge but i dont understand what orientation the circles (created by the closed rings) have. i will look into the cw etc.
s1 = Sphere((0.0,0.0),1.0)
s2 =s1 |> Scale(2)
r1 = discretize(s1)
viz(r1, showsegments = true)
r2 = discretize(s2)
p = PolyArea(r2,reverse(r1))
m = discretize(p) |> viz
ERROR: MethodError: no method matching reverse(::SimpleMesh{𝔼{…}, CoordRefSystems.Cartesian2D{…}, Vector{…}, GridTopology{…}})
The function `reverse` exists, but no method is defined for this combination of argument types.
Closest candidates are:
reverse(::PlotUtils.CategoricalColorGradient)
@ PlotUtils ~/.julia/packages/PlotUtils/dVEMd/src/colorschemes.jl:152
reverse(::PlotUtils.ColorPalette)
@ PlotUtils ~/.julia/packages/PlotUtils/dVEMd/src/colorschemes.jl:265
reverse(::PlotUtils.ContinuousColorGradient)
@ PlotUtils ~/.julia/packages/PlotUtils/dVEMd/src/colorschemes.jl:76
...
Stacktrace:
[1] top-level scope
@ Untitled-2:35
Some type information was truncated. Use `show(err)` to see complete types.
I did get some further by reversing the direction of the Y vector I put into tuples. Any idea of how I can refine the mesh further and bring more balanced triangles in?
# reference sphere
s = Sphere((0, 0), 1)
# parameter range
ts = 0.01:0.01:0.99
# outer ring
o = Ring([s(t) for t in ts])
# inner ring
i = Ring([s(t) |> Scale(0.5) for t in reverse(ts)])
# polygonal area
p = PolyArea(o, i)
m = discretize(p, DelaunayTriangulation())
viz(m, showsegments=true)
These triangles are not good for the FEM. Recent releases of DelaunayTriangulation.jl allow discretization with outer and inner boundaries, so we should probably update our code to produce better triangles by default.
What you can do currently is refine the mesh with TriRefinement and a predicate function, or use DelaunayTriangulation.jl directly.
@Br0BOLE the Delaunay triangulation takes a set of points as input and a set of indices of points marking the boundary. By default we only use the boundary points of polygons, leading to the triangles illustrated above.
If you need triangles with better angles, you can insert more points in the list. Below is a helper function that samples points inside the polygonal area according to a minimum distance, and then calls the underlying Delaunay routine to create the mesh:
using Meshes
using Unitful
import DelaunayTriangulation as DT
# alpha is the minimum distance between points
function mydiscretize(poly::Polygon, α=0.05)
verts = map(rings(poly)) do ring
collect(eachvertex(ring))
end
offset = 0
bverts = [[Int[]] for _ in 1:length(verts)]
for i in 1:length(verts)
nverts = length(verts[i])
bverts[i][1] = [(offset + 1):(offset + nverts); (offset + 1)]
offset += nverts
end
points1 = reduce(vcat, verts)
points2 = sample(poly, MinDistanceSampling(α)) |> collect
points = [points1; points2]
coords = map(p -> ustrip.(to(p)), points)
triang = DT.triangulate(coords, boundary_nodes=bverts)
connec = connect.(DT.each_solid_triangle(triang))
SimpleMesh(points, connec)
end
m = mydiscretize(p, 0.05)
viz(m, showsegments=true)
This function is a workaround. We will try to find the time to incorporate this function into our Meshes.jl machinery to faciliate the lives of end-users.
The downside of the above implementation is that it relies on another discretization of the geometry to sample points in the interior. We can do better with primitive geometries like balls, sphere, etc.
I ran the example with the same polygon p of my previous comment above.
