# FE Discretisation for embedded cables

**URL:** <https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025>\
**Category:** Geo\
**Tags:** meshes\
**Created:** [March 16, 2025, 9:14pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025 "2025-03-16T21:14:07Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 16, 2025, 9:14pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/1 "2025-03-16T21:14:07Z")

</div>

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?

```julia
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

```

Many thanks for any suggestions

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 16, 2025, 9:21pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/2 "2025-03-16T21:21:59Z")

</div>

Meshes.jl maintainer here. What is the issue you’re having ?

Consider `Circle`, `Disk`, `Sphere` and `Ball` as related primitive geometries.

You can discretize these geometries to get rings, meshes etc.

Notice that rings of a polygon must have a certain orientation. The outer ring must be CCW and the inner rings must be CW.

You can fix that in general by sending the `PolyArea` to the `Repair(11)` transform.

---

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 16, 2025, 9:33pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/3 "2025-03-16T21:33:11Z")

</div>

the mesh created by the example above doesnt seem to fully respect the boundary (circular hole in the middle).

i want to assign different thermal parameters inside and outside the hole created. does that make sense?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 16, 2025, 9:33pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/4 "2025-03-16T21:33:56Z")

</div>

Can you please make sure the orientation of rings is correct as explained above?

---

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 16, 2025, 9:37pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/5 "2025-03-16T21:37:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 16, 2025, 9:42pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/6 "2025-03-16T21:42:14Z")

</div>

No worries. I am typing from my phone, but try the following:

```julia
s1 = Sphere((0,0), 1)
s2 = s1 |> Scale(2)

r1 = discretize(s1)
r2 = discretize(s2)

p = PolyArea(r2, reverse(r1))

m = discretize(p)

```

---

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 16, 2025, 10:01pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/7 "2025-03-16T22:01:25Z")

</div>

```julia

```

```julia
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?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 16, 2025, 10:39pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/8 "2025-03-16T22:39:54Z")

</div>

I dont recall now but maybe we changed the output type of discretize on 2D sphere.

You can check the official docs of the package to learn about the different discretization methods available. It should contain examples.

The reverse of the inner ring is all you need to create a valid PolyArea.

The discretization methods will take the PolyArea as input and will produce the mesh with the properties you are after.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 17, 2025, 11:17am UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/9 "2025-03-17T11:17:36Z")

</div>

@Br0BOLE below is a commented example:

```julia
# 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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/a/0a0b1ae6517ab8b6da6969ed7f0ad8fb0ff40409.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 17, 2025, 11:39am UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/10 "2025-03-17T11:39:21Z")

</div>

You can also just use an external mesh generator like Gmsh, and load the resulting mesh into Julia. Or use it directly from Julia, via Gmsh.jl.

---

<div class="post-metadata">

**Author:** ![j-fu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j-fu/32/11373_2.png) [@j-fu](https://discourse.julialang.org/u/j-fu)\
**Post date:** [March 17, 2025, 11:46am UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/11 "2025-03-17T11:46:15Z")

</div>

FYI, mesh generators TetGen and Triangle are available via TetGen.jl and Triangulate.jl, respectively.

We also maintain ExtendableGrids.jl and SimplexGridFactory.jl to ease their use.

---

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 17, 2025, 11:53am UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/12 "2025-03-17T11:53:04Z")

</div>

thanks very much all. i will contimue to play with meshes and the other suggested package and get back with any conclusions Juliohm.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [March 18, 2025, 12:31pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/13 "2025-03-18T12:31:52Z")

</div>

@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:

```julia
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)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce468bc69617760dcda5571e2be884a59a6538c8.jpeg)

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.

---

<div class="post-metadata">

**Author:** ![Br0BOLE](https://avatars.discourse-cdn.com/v4/letter/b/b19c9b/32.png) [@Br0BOLE](https://discourse.julialang.org/u/Br0BOLE)\
**Post date:** [March 18, 2025, 10:12pm UTC](https://discourse.julialang.org/t/fe-discretisation-for-embedded-cables/127025/14 "2025-03-18T22:12:38Z")

</div>

thanks Juliohm