Creating an ellipse mesh produces spurious points

[griff10000]

Creating an ellipse generates spurious points on the periphery. Windows 10, python 3.8.3, Anaconda.

import dmsh
geo = dmsh.Ellipse([0.0, 0.0], 2.0, 1.0)

edge_size=0.2

X, cells = dmsh.generate(geo, edge_size, tol=1.0e-10) # changing tol does not help

dmsh.helpers.show(X, cells, geo)

[oddroj]

This was a problem that I encountered semi-recently. I found a workaround that produced a mesh that still was not very good, but was better for my application.

This mesh, above, is created using a constant angular spacing of nodes initially, essentially you are meshing in polar coordinates. This is why it bunches up in places. If you can eliminate that problem then creating a more usable mesh is easier. To do this you need to use the pseudo-angle from the elliptical coordinate system.

Then you mesh initially over the rectangular in elliptical coordinates domain, then you can convert the mesh to the cartesian coordinate system.

import numpy as np
import optimesh
import dmsh
import matplotlib.pyplot as plt

def EllipToCart(Nu,Zeta,A,C):
    #Nu is the pseudo  angular elliptical coordinate
    # Zeta is the pseudo radial elliptical coordinate
    # A is the elliptical radius in the Y direction
    # C is the elliptical radius in the X direction
    if A/C<1:
        b=np.sqrt(np.abs(A**2 - C**2))
        x=b*np.cosh(Zeta)*np.cos(Nu)
        y=b*np.sinh(Zeta)*np.sin(Nu)

        Jacobian= 1/2 * b**2 *( np.cosh(2*Zeta)-np.cos(2*Nu))

    elif A/C>1:
        b=np.sqrt(np.abs(A**2 - C**2))
        x=b*np.sinh(Zeta)*np.cos(Nu)
        y=b*np.cosh(Zeta)*np.sin(Nu)
        Jacobian = 1/2 * b**2 *( np.cosh(2*Zeta)+np.cos(2*Nu))
    elif A/C==1:
        x=Zeta*np.cos(Nu)
        y=Zeta*np.sin(Nu)
        Jacobian=Zeta

    return x,y, Jacobian

def EllipMesh(res,A,C,elliptical_coords):
     b=np.sqrt(np.abs(A**2 - C**2))
     if A==C:
          Zmax=A
     else:
          Zmax=np.log((A+C)/b)
     Nmax=2*np.pi

     geo=dmsh.Rectangle(0,Nmax,0,Zmax)
     points_ellipse, cells_ellipse= dmsh.generate(geo, min(Nmax,Zmax)/res)

     if not elliptical_coords:
         geo=dmsh.Ellipse([0,0],C,A)

         x,y,j = EllipToCart(points_ellipse[:,0], points_ellipse[:,1],A,C)
         points=np.array(np.array([x,y]).T)

         points, cells = optimesh.optimize_points_cells(
             points, cells_ellipse, "CVT (full)", 1.0e-10, 100)
         return points, cells, geo

     else:
             points, cells = optimesh.optimize_points_cells(
                 points_ellipse, cells_ellipse, "CVT (full)", 1.0e-10, 100)
             return points, cells, geo

plt.figure()

A=4
C=10

P,C,G = EllipMesh(10,A,C,0)
dmsh.helpers.show(P,C,G)

Gives.

Note the displaced spurious points shown up close here.

[griff10000]

Many thanks. Whilst not perfect, this is a good improvement.

[griff10000]

I have thought a bit more about this problem and come up with a solution that generates almost equidistance points around an ellipse, and then uses the polygon mesh feature in dmsh to create a mesh over the ellipse. I adapted a piece of python code from the web that generates the points on the ellipse boundary. The code is as follows:

# Generate equidistance (almost) points around the periphery of an ellipse
# Ref: https://stackoverflow.com/questions/44334655/how-to-divide-circumference-of-an-ellipse-equally?noredirect=1&lq=1

# Note: A more exact procedure is possible using 'elliptic integrals'. Finding
#       the true value for the perimiter is straight forward using python
#       function ellipe and is included as a check. However, a full calculation
#       requires obtaining inverse elliptic integrals and I have not yet found
#       a way yet. So have abandoned this for now., see:
# Ref: https://math.stackexchange.com/questions/172766/calculating-equidistant-points-around-an-ellipse-arc

import numpy as np
from scipy.special import ellipe # complete elliptic integral of the second type
import matplotlib.pyplot as plt

def distance(x1,y1,x2,y2):
    return np.sqrt((x2-x1)**2 + (y2-y1)**2)

a = 10 #5
b = 4 #3
m = 1-b**2/a**2    # eccentricity of ellipse
da = 4*a*ellipe(m) # analytical value for ellipse perimiter

x0 = a
y0 = 0
angle = 0
d = 0
del_a = 2*np.pi/10000
while(angle<=2*np.pi):
  x = a * np.cos(angle)
  y = b * np.sin(angle)
  d += distance(x0,y0,x,y)
  x0 = x
  y0 = y
  angle += del_a
print( "Circumference of ellipse: numeric = %f, analytic = %f" %(d,da))
N = 17
h = d/N
angle = 0
x0 = a
y0 = 0
angle0 = 0
xx = np.zeros(N)
yy = np.zeros(N)
for i in range(N):
  dist = 0
  while(dist<h):
    x = a * np.cos(angle)
    y = b * np.sin(angle)
    dist += distance(x0,y0,x,y)
    x0 = x
    y0 = y
    angle += del_a
  xx[i] = x
  yy[i] = y

# generate analytical points on ellipse 
N1 = 50
aa = np.linspace(0,2*np.pi,N1)
xa = a*np.cos(aa)
ya = b*np.sin(aa)

# plt.figure()
# plt.subplot(211)
# plt.plot(xa,ya,"b-")
# plt.plot(xx,yy,'r.')

# =================
# Use dmsh to generate mesh from polygon

import dmsh
import optimesh

# ellipse data as a polygon with N1 points
dat = [[xa[0],ya[0]]]
for i in range(1,N1-1) :
  dat.append([xa[i],ya[i]])

geo = dmsh.Polygon(dat)

#plt.subplot(212)

X, cells = dmsh.generate(geo, 0.5)

# optionally optimize the mesh
import optimesh
X, cells = optimesh.optimize_points_cells(X, cells, "CVT (full)", 1.0e-10, 100)
# visualize the mesh
dmsh.helpers.show(X, cells, geo, full_screen=False)

I have not tried to optimise the code. It generates: and

[griff10000]

I have now managed to come up with a solution that generates equidistance points around an ellipse using elliptic functions, which then uses the polygon mesh feature in dmsh to create a mesh over the ellipse. The code below runs fast on my windows PC:

# Calculate equidistance points around the periphery of an ellipse and then use
# the 'dmsh' polygon feature to generate a mesh over the ellipse.

# Note: An exact procedure is possible using 'elliptic integrals'. Finding
#       the true value for the perimiter is straight forward using python
#       function 'ellipe()'. However, a full calculation requires obtaining
#       inverse elliptic integrals. These are not directly available as python
#       callable functions, so I have used python functions 'fsolve()' and
#       'ellipeinc()' to obtain the required values. See also:
# Ref: https://math.stackexchange.com/questions/172766/calculating-equidistant-points-around-an-ellipse-arc

import numpy as np 
import matplotlib.pyplot as plt
import scipy.optimize
# load functions for complete and incomplete elliptic integrals
# of the second type
from scipy.special import ellipe, ellipeinc 

# Set up parameters for ellipse
a = 5
b = 3
m = 1-b**2/a**2   # eccentricity of ellipse
d = 4*a*ellipe(m) # analytical total length of ellipse perimiter

N = 51            # number of points on perimiter

di = np.linspace(0, d, N)        # Equally spaced points around perimiter
dphi = np.linspace(0,2*np.pi, N) # starting points for fsolve function

# function used with fsolve
def fun(phi, data) :
  d, m = data
  res = d - a*ellipeinc(phi, m)
  return res

# Solve for angles defining point around ellipse
mi = m*np.ones(N) # vectorise 'm'
phi= scipy.optimize.fsolve(fun,dphi, xtol=1e-12, args=[di,mi])

# calculate (x,y) values for points around perimiter
x = a * np.sin(phi)
y = b * np.cos(phi)

# Plot ellipse with points superimposed
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)
ax.plot(x,y,'b-')
ax.plot(x,y,'r.')
ax.set_aspect('equal')
ax.grid()

# =================
# Use dmsh to generate mesh from polygon

import dmsh
import optimesh

# ellipse data as a polygon with N1 points
dat = [[x[0],y[0]]]
for i in range(1,N-1) :
  dat.append([x[i],y[i]])

geo = dmsh.Polygon(dat)

plt.figure(figsize=(10, 8))

X, cells = dmsh.generate(geo, 0.25)

# optionally optimize the mesh
X, cells = optimesh.optimize_points_cells(X, cells, "CVT (full)", 1.0e-10, 100)
# visualize the mesh
dmsh.helpers.show(X, cells, geo, full_screen=False)

the code generates: and