Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Mapping a disc to a square using FG method

Tags:

python

math

mesh

Context:

I am working with meshes in computational fluid dynamics. I would like to generate a structured mesh around a circle. My plan is to generate a polar mesh (Left mesh in Figure below) then use FG formulas to get the final mesh (Right mesh on Figure below). I am using FG method from Page 4 of this article to map a disc to square . Unfortunately the article doesn't mention how to handle singularities in the formulas. Here are the expressions:

x = sgn(uv)/(v*sqrt 2)*sqrt(u**2+v**2-sqrt((u**2+v**2)(u**2+v**2-4u**2v**2)))
y = sgn(uv)/(u*sqrt 2)*sqrt(u**2+v**2-sqrt((u**2+v**2)(u**2+v**2-4u**2v**2)))

enter image description here

Before programming this, I am struggling with some problems with these formaulas

Questions

  • Why do these formulae map the following points: (1,0), (0,1), (-1,0), (-1,-1), and (0,0) to the point (0,0) ?
  • How I am supposed to get the intermediate shape between a circle and a square as shown in the figure below.

  • Is it possible to provide an algorithm to get the right mesh from the left one?

enter image description here

Here is my attempt:

"""Map a circular computational domain with structured mesh around a circle (circular cylinder in 3D) to
Rectangular domain"""

import numpy as np
from numpy import sqrt, sign, pi, cos, sin
import matplotlib.pyplot as plt

def FGsquircle(u, v):
    SMALL = 1e-15
    t0 = u**2+v**2
    t1 = (u**2+v**2)*(u**2+v**2-4*u**2*v**2)
    t2 = u**2+v**2
    t3 = (u**2+v**2)*(u**2+v**2-4*u**2*v**2)

    x = sign(u*v)/(v*sqrt(2.0)+SMALL)*sqrt(t0-sqrt(t1))
    y = sign(u*v)/(u*sqrt(2.0)+SMALL)*sqrt(t2-sqrt(t3))
    return x, y

R0 = 1.0 # radius of the disc
RMAX = 5.0 # the radius of the outer circle in the domain
NT = 360 # num of division in the theta direction
NR = 10 # num of radial divisions
r = [R0+(RMAX-R0)/NR*k for k in range(NR)] # the radii of circles
theta = np.array([2*pi/NT*k for k in range(NT+1)])

u = [r[k]*cos(theta) for k in range(NR)]
v = [r[k]*sin(theta) for k in range(NR)]

u = np.array(u)
v = np.array(v)

x, y = FGsquircle(u, v)

I got the following error:

utils.py:21: RuntimeWarning: invalid value encountered in sqrt
  x = sign(u*v)/(v*sqrt(2.0)+SMALL)*sqrt(t0-sqrt(t1))
utils.py:22: RuntimeWarning: invalid value encountered in sqrt
  y = sign(u*v)/(u*sqrt(2.0)+SMALL)*sqrt(t2-sqrt(t3))

I appreciate any help.


1 Answers

Why do these formulae map the following points: (1,0), (0,1), (-1,0), (-1,-1), and (0,0) to the point (0,0) ?

In both expressions, you divide by u or v, so when one of the two is 0, the expression becomes undefined (not zero).

How I am supposed to get the intermediate shape between a circle and a square as shown in the figure below.

Just map a circle of radius smaller than 1.

Is it possible to provide an algorithm to get the right mesh from the left one?

You'll just need to transform the mesh points. The cells can remain the same.


Example code:

from numpy import sqrt, sign
import numpy
import matplotlib.pyplot as plt


def f(x):
    u, v = x
    alpha = sqrt(
        u ** 2
        + v ** 2
        - sqrt((u ** 2 + v ** 2) * (u ** 2 + v ** 2 - 4 * u ** 2 * v ** 2))
    )
    return numpy.array(
        [sign(u * v) / (v * sqrt(2)) * alpha, sign(u * v) / (u * sqrt(2)) * alpha]
    )


for r in numpy.linspace(0.1, 1.0, 10):
    theta = numpy.linspace(0.0, 2 * numpy.pi, 1000, endpoint=True)
    uv = r * numpy.array([numpy.cos(theta), numpy.sin(theta)])
    xy = f(uv)

    plt.plot(xy[0], xy[1], "-")

plt.gca().set_aspect("equal")
plt.show()

enter image description here

like image 156
Nico Schlömer Avatar answered Nov 30 '25 11:11

Nico Schlömer



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!