3 votes

Cartographie d'un disque en un carré à l'aide de la méthode FG

Contexte :

Je travaille avec des maillages dans le domaine de la dynamique des fluides. J'aimerais générer un maillage structuré autour d'un cercle. Mon plan consiste à générer un maillage polaire (maillage de gauche dans la figure ci-dessous), puis à utiliser les formules FG pour obtenir le maillage final (maillage de droite dans la figure ci-dessous). J'utilise la méthode FG de Page 4 of cet article pour faire correspondre un disque à un carré . Malheureusement, l'article ne mentionne pas comment traiter les singularités dans les formules. Voici les 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

Avant de programmer cela, j'ai quelques problèmes avec les formules suivantes

Questions

  • Pourquoi ces formules représentent-elles les points suivants : (1,0) , (0,1) , (-1,0) , (-1,-1) et (0,0) jusqu'au bout (0,0) ?

  • Comment suis-je censé obtenir la forme intermédiaire entre un cercle et un carré comme indiqué dans la figure ci-dessous.

  • Est-il possible de fournir un algorithme permettant d'obtenir le bon maillage à partir du maillage de gauche ?

enter image description here

Voici ma tentative :

"""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)

J'ai obtenu l'erreur suivante :

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

J'apprécie toute aide.

2voto

Nico Points 893

Pourquoi ces formules font-elles correspondre les points suivants : (1,0), (0,1), (-1,0), (-1,-1) et (0,0) au point (0,0) ?

Dans les deux expressions, vous divisez par u o v Ainsi, lorsque l'un des deux est égal à 0, l'expression devient indéfini (pas zéro).

Comment suis-je censé obtenir la forme intermédiaire entre un cercle et un carré comme indiqué dans la figure ci-dessous.

Il suffit de tracer un cercle de rayon inférieur à 1.

Est-il possible de fournir un algorithme permettant d'obtenir le bon maillage à partir du maillage de gauche ?

Il suffit de transformer les points du maillage. Les cellules peuvent rester inchangées.


Exemple de 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

Prograide.com

Prograide est une communauté de développeurs qui cherche à élargir la connaissance de la programmation au-delà de l'anglais.
Pour cela nous avons les plus grands doutes résolus en français et vous pouvez aussi poser vos propres questions ou résoudre celles des autres.

Powered by:

X