How to find the conformal mapping between two regions

Question

As the Riemann Mapping Theorem, if both regions are simply connected regions, there must be a unique one-conformal mapping such that the points between the two regions correspond one-to-one. Can we find a conformal mapping from the red region to the blue region using MMA?

Show[RegionPlot[x^2 + y^3 < 2 && y > x^2, {x, -2, 3}, {y, -1, 4}, 
  PlotStyle -> Red], 
 Graphics[{Blue, Triangle[{{2, 1}, {2, 3}, {0, 4}}]}]]

ps: I can accept a polygon that discretes the arc into many sides.

Isn't this answer mathematica.stackexchange.com/a/178333/9469 completely solves your problem? — yarchik, Commented Jan 25 at 22:13
Sorry, unfortunately, I have no time for a detailed answer. Fortunately, there a people who have worked much more on this topic. Please have a look into this article by Mark Gillespie, Boris Springborn, and Keenan Crane: dl.acm.org/doi/10.1145/3450626.3459763 and into the references therein. (There is also a video of Mark's talk.) I know, this is not the anwer you were looking for. But it is the best that I can give to you at the moment. — Henrik Schumacher, Commented Jan 26 at 12:12
Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematica Meta, or in Mathematica Chat. Comments continuing discussion may be removed. — Kuba, Commented Jan 31 at 12:41

Alex Trounev · Accepted Answer · 2024-02-07 14:55:29Z

We can use PI algorithm described in the paper Numerical Conformal Mapping with Rational Functions.The function for conformal mapping can be computed as follows

Needs["NDSolve`FEM`"];
reg1 = Triangle[{{2, 1}, {2, 3}, {0, 4}}]; reg2 = 
 ImplicitRegion[x^2 + y^3 < 2 && y > x^2, {x, y}];

f[reg_, order_] := 
  Module[{Z, A, B, eq, g, vec, mat, sol, R1, z1, bndedges, 
    bndvertices, bpts, b, b0, aa, bb}, z1 = N@RegionCentroid[reg]; 
   R1 = ToElementMesh[reg, "MaxBoundaryCellMeasure" -> .01, 
     "MaxCellMeasure" -> .001]; 
   bndedges = R1["BoundaryElements"][[1, 1]];
   bndvertices = Sort@DeleteDuplicates[Flatten[bndedges]];
   bpts = R1["Coordinates"][[bndvertices]];
   b = ConstantArray[0., Length[R1["Coordinates"]]];
   b0 = b[[
      bndvertices]] = -0.5 Log[
       Total[(bpts - ConstantArray[z1, Length[bpts]])^2, {2}]];
   
   With[{nn = order}, 
    Z = {x, y} |-> Table[ComplexExpand[(x + I y)^n], {n, 0, nn}];
    A = Array[aa, {nn + 1}]; B = Array[bb, {nn + 1}];
    g = {x, y} |-> ComplexExpand[(A + I B) . Z[x, y]]];
   
   eq = -b0 + (g[#[[1]] - z1[[1]], #[[2]] - z1[[2]]] & /@ bpts) /. 
     I -> 0;
   
   {vec, mat} = CoefficientArrays[Join[{bb[1]}, eq], Join[A, B]];
   
   sol = LeastSquares[mat, -vec];
   With[{nn = order}, 
    Do[aa[i] = sol[[i]]; 
     bb[i] = sol[[nn + 1 + i]];, {i, nn + 1}]]; {{x, y} |-> 
     ReIm[Evaluate[
       ComplexExpand[((x + I y) - (z1[[1]] + I  z1[[2]]))  Exp[
          g[x - z1[[1]], y - z1[[2]]]]]]], {x, y} |-> 
     ComplexExpand[(A + I B) . Z[x, y]], z1, R1}];

Example of usage

{f1, g1, z1, mesh1} = f[reg1, 32];

{Plot3D[Evaluate[g1[x - z1[[1]], y - z1[[2]]]] // Re, 
  Element[{x, y}, reg1], ColorFunction -> Hue], 
 Plot3D[Evaluate[g1[x - z1[[1]], y - z1[[2]]]] // Im, 
  Element[{x, y}, reg1], ColorFunction -> Hue]}

{f2, g2, z2, mesh2} = f[reg2, 32];

{Plot3D[Evaluate[g2[x - z2[[1]], y - z2[[2]]]] // Re, 
  Element[{x, y}, reg2], ColorFunction -> Hue], 
 Plot3D[Evaluate[g2[x - z2[[1]], y - z2[[2]]]] // Im, 
  Element[{x, y}, reg2], ColorFunction -> Hue]}

Mesh lines on reg2

rho2 = ContourPlot[Abs[f2[x, y]], Element[{x, y}, reg2], 
   Contours -> Range[.1, .9, .1], 
   ColorFunction -> Function[{x, y}, Hue[1/3, 1, 1, .5]], 
   AspectRatio -> Automatic];

arg2 = ContourPlot[Arg[f2[x, y]], Element[{x, y}, reg2], 
   Contours -> Pi  Range[-.8, .8, .1], 
   ColorFunction -> Function[{x, y}, Hue[1/3, 1, 1, .25]], 
   AspectRatio -> Automatic, ExclusionsStyle -> Black];

Show[rho2, arg2]

Mesh lines on reg1

rho1 = ContourPlot[Abs[f1[x, y]], Element[{x, y}, reg1], 
  Contours -> Range[.1, .9, .1], 
  ColorFunction -> Function[{x, y}, Hue[1/3, 1, 1, .5]], 
  AspectRatio -> Automatic];

arg1 = ContourPlot[Arg[f1[x, y]], Element[{x, y}, reg1], 
  ColorFunction -> Function[{x, y}, Hue[1/3, 1, 1, .25]], 
  AspectRatio -> Automatic, Contours -> Pi Range[-.8, .8, .1], 
  PlotPoints -> 50];
Show[rho1, arg1]

First attempt. We can try approach that recommended by yarchik and implemented by
Henrik Schumacher here. We have

Needs["NDSolve`FEM`"];
reg1 = Triangle[{{2, 1}, {2, 3}, {0, 4}}]; reg2 = 
 ImplicitRegion[x^2 + y^3 < 2 && y > x^2, {x, y}];
f[reg_] := 
 Module[{R, p, z0, ufun, vfun}, 
  R = ToElementMesh[reg, "MeshOrder" -> 1, 
    "MaxBoundaryCellMeasure" -> .01, "MaxCellMeasure" -> .001];
  z0 = N@RegionCentroid[reg];
  (*Initialization of Finite Element Method*)
  vd = NDSolve`VariableData[{"DependentVariables", 
      "Space"} -> {{u}, {x, y}}];
  sd = NDSolve`SolutionData[{"Space"} -> {R}];
  cdata = 
   InitializePDECoefficients[vd, sd, 
    "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}, 
    "MassCoefficients" -> {{1}}];
  bcdata = 
   InitializeBoundaryConditions[vd, 
    sd, {DirichletCondition[u[x, y] == 0., True]}];
  mdata = InitializePDEMethodData[vd, sd];
  
  (*Discretization*)
  dpde = DiscretizePDE[cdata, mdata, sd];
  dbc = DiscretizeBoundaryConditions[bcdata, mdata, sd];
  {load, stiffness, damping, mass} = dpde["All"];
  DeployBoundaryConditions[{load, stiffness}, dbc];
  
  (*Preparation of Dirichlet boundary conditions for u*)
  
  bndedges = R["BoundaryElements"][[1, 1]];
  bndvertices = Sort@DeleteDuplicates[Flatten[bndedges]];
  bpts = R["Coordinates"][[bndvertices]];
  b = ConstantArray[0., Length[R["Coordinates"]]];
  b[[bndvertices]] = -0.5 Log[
     Total[(bpts - ConstantArray[z0, Length[bpts]])^2, {2}]];
  
  (*Solving the system and creating an interpolating function*)
  
  solver = LinearSolve[stiffness, Method -> "Pardiso"];
  uvals = solver[b];
  ufun = ElementMeshInterpolation[{R}, uvals];
  (*Preparation of Neumann boundary conditions for v*)
  gradu = {x, y} |-> Evaluate[D[ufun[x, y], {{x, y}, 1}]];
  J = N@RotationMatrix[Pi/2];
  p = R["Coordinates"];
  {i, j} = Transpose[R["BoundaryElements"][[1, 1]]];
  normalprojections = 
   MapThreadDot[
    R["BoundaryNormals"][[
     1]], (gradu @@@ (0.5 (p[[i]] + p[[j]]))) . (-J)];
  boundaryedgelengts = Sqrt[Total[(p[[i]] - p[[j]])^2, {2}]];
  {\[Alpha], \[Beta]} = Transpose[bndedges];
  vertexbndedgeconnectivity = 
   SparseArray[
    Transpose[{Join[\[Alpha], \[Beta]], 
       Join[Range[Length[\[Alpha]]], Range[Length[\[Beta]]]]}] -> 
     1, {Length[p], Length[bndedges]}];
  
  (*Solving the system and creating an interpolating function*)
  
  b = vertexbndedgeconnectivity . (normalprojections \
boundaryedgelengts);
  vvals = solver[b];
  vfun = ElementMeshInterpolation[{R}, vvals];
  {x, y} |-> 
   Evaluate[
    ComplexExpand[
     ReIm[((x + I y) - (z0[[1]] + I z0[[2]])) Exp[
        ufun[x, y]] (Cos[vfun[x, y]] + I  Sin[vfun[x, y]])]]]];

We use function f at reg1, reg2 as follows

f1 = f[reg1];

With[{R = 
   ToElementMesh[reg1, "MeshOrder" -> 1, 
    "MaxBoundaryCellMeasure" -> .01, "MaxCellMeasure" -> .001]}, 
 p = R["Coordinates"];
 tex = ColorData["BrightBands", "Image"];
 texcoords = 
  Transpose[{Total[(f1 @@@ p)^2, {2}], ConstantArray[0.5, Length[p]]}];
 g1 = {Graphics[{Texture[tex], 
     ElementMeshToGraphicsComplex[R, 
      VertexTextureCoordinates -> texcoords]}, PlotRange -> All], 
   Graphics[{Texture[tex], 
     ElementMeshToGraphicsComplex[R, 
      VertexTextureCoordinates -> texcoords, 
      "CoordinateConversion" -> (f1 @@@ # &)]}]}]

f2 = f[reg2];

With[{R = 
   ToElementMesh[reg2, "MeshOrder" -> 1, 
    "MaxBoundaryCellMeasure" -> .01, "MaxCellMeasure" -> .001]}, 
 p = R["Coordinates"];
 tex = ColorData["BrightBands", "Image"];
 texcoords = 
  Transpose[{Total[(f2 @@@ p)^2, {2}], ConstantArray[0.5, Length[p]]}];
 g1 = {Graphics[{Texture[tex], 
     ElementMeshToGraphicsComplex[R, 
      VertexTextureCoordinates -> texcoords]}, PlotRange -> All], 
   Graphics[{Texture[tex], 
     ElementMeshToGraphicsComplex[R, 
      VertexTextureCoordinates -> texcoords, 
      "CoordinateConversion" -> (f2 @@@ # &)]}]}]

Nothing new above, but now we can map reg2 to reg1 using FindRoot

z0 = N@RegionCentroid[reg2];

ps = Sort[p, Norm[# - z0] &];

ds = f2[#[[1]], #[[2]]] & /@ ps;

z1 = N@RegionCentroid[reg1];
ps1 = Table[{x, y} /. 
    Quiet@FindRoot[
      f1[x, y] == ds[[i]], {x, z1[[1]]}, {y, z1[[2]]}], {i, 
    Length[ds]}];

Visualization

kmax = First[Last[Length[ps] // FactorInteger]];

nmax = Length[ps]/kmax;

{ListPlot[
  Table[Take[ps, {kmax  i + 1, kmax  (i + 1)}], {i, 0, nmax - 1}], 
  PlotStyle -> cl], 
 ListPlot[
  Table[Take[ps1, {kmax  i + 1, kmax  (i + 1)}], {i, 0, nmax - 1}], 
  AspectRatio -> 1, PlotStyle -> cl, PlotLegends -> Automatic]}

We also can check how parts of reg2 map at reg

cl = {Red, Orange, Green, Blue, Black, Gray};

Table[{ListPlot[Take[ps, {kmax  i + 1, kmax  (i + 1)}], 
   PlotStyle -> cl[[i + 1]]], 
  ListPlot[Take[ps1, {kmax  i + 1, kmax  (i + 1)}], AspectRatio -> 1, 
   PlotStyle -> cl[[i + 1]]]}, {i, 0, nmax - 1}]

Update 1. Given the pictures that the user cvgmt provided us (thanks to him), we need to check how Schumacher's code calculates the functions ufun, vfun. For this we used NDSolve and our code with PI algorithm implemented in it - see paper Numerical Conformal Mapping with Rational Functions. NDSolve code is very simple

Needs["NDSolve`FEM`"];
reg1 = Triangle[{{2, 1}, {2, 3}, {0, 4}}]; reg2 = 
 ImplicitRegion[x^2 + y^3 < 2 && y > x^2, {x, y}];

z1 = N@RegionCentroid[reg1]; R1 = 
 ToElementMesh[reg1, "MaxBoundaryCellMeasure" -> .01, 
  "MaxCellMeasure" -> .001];
U1 = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 0, 
    DirichletCondition[
     u[x, y] == -.5 Log[(x - z1[[1]])^2 + (y - z1[[2]])^2], True]}, u,
    Element[{x, y}, R1]];
Plot3D[U1[x, y], Element[{x, y}, reg1], ColorFunction -> Hue]

The PI algorithm is based on the series expansion of a function $g(z)=u(z)+i v(z)$, where $\nabla^2 u=0$ with the Dirichlet condition $u=-\ln|z-z1|$. To compute $u,v$ we use mesh R1 generated above:

bndedges = R1["BoundaryElements"][[1, 1]];
bndvertices = Sort@DeleteDuplicates[Flatten[bndedges]];
bpts = R1["Coordinates"][[bndvertices]];
b = ConstantArray[0., Length[R1["Coordinates"]]];
b0 = b[[bndvertices]] = -0.5 Log[
     Total[(bpts - ConstantArray[z1, Length[bpts]])^2, {2}]];

With[{nn = 32}, 
  Z = {x, y} |-> Table[ComplexExpand[(x + I y)^n], {n, 0, nn}]; 
  A = Array[aa, {nn + 1}]; B = Array[bb, {nn + 1}]; 
  g = {x, y} |-> ComplexExpand[(A + I B) . Z[x, y]]];

eq = -b0 + (g[#[[1]] - z1[[1]], #[[2]] - z1[[2]]] & /@ bpts) /. I -> 0;

{vec, mat} = CoefficientArrays[Join[{bb[1]}, eq], Join[A, B]];

sol = LeastSquares[mat, -vec];
With[{nn = 32}, 
  Do[aa[i] = sol[[i]]; bb[i] = sol[[nn + 1 + i]];, {i, nn + 1}]];

Now we can compare U1 computed with NDSolve, ufun computed with Henrik Schumacher code (function f above), and function u=Re[g[z-z1]]

{Plot3D[U1[x, y], Element[{x, y}, reg1], ColorFunction -> Hue], 
 Plot3D[ufun[x, y], Element[{x, y}, reg1], ColorFunction -> Hue], 
 Plot3D[Evaluate[g[x - z1[[1]], y - z1[[2]]]] // Re, 
  Element[{x, y}, reg1], ColorFunction -> Hue]}

Fortunately all 3 functions look same

But unfortunately functions vfun and Im[g[z-z1]] are differ

{Plot3D[vfun[x, y], Element[{x, y}, reg1], ColorFunction -> Hue], 
 Plot3D[Evaluate[g[x - z1[[1]], y - z1[[2]]]] // Im, 
  Element[{x, y}, reg1], ColorFunction -> Hue]}

Using g[z] computed above we can define conformal mapping reg1 on the unit disk as follows

f1p = {x, y} |-> 
   ReIm[Evaluate[
     ComplexExpand[((x + I y) - (z1[[1]] + I  z1[[2]]))  Exp[
        g[x - z1[[1]], y - z1[[2]]]] ]]];

Finally we test f1p using code provided by cvgmt

disktoreg1 = 
  ParametricPlot[{x, y}, {x, y} \[Element] reg1, 
   MeshFunctions -> {Function[{x, y}, 
      Norm@{Indexed[f1p[x, y], 1], Indexed[f1p[x, y], 2]}], 
     Function[{x, y}, 
      ArcTan @@ {Indexed[f1p[x, y], 1], Indexed[f1p[x, y], 2]}]}, 
   Frame -> False, Axes -> False, Mesh -> {30, 30}];
reg1todisk = 
  ParametricPlot[f1p[x, y], {x, y} \[Element] reg1, 
   MeshFunctions -> {Function[{x, y}, Norm@{x, y}], 
     Function[{x, y}, ArcTan @@ {x, y}]}, Frame -> False, 
   Axes -> False, Mesh -> {30, 30}];
GraphicsRow[{disktoreg1, reg1todisk}]

Now it looks like conformal map since all grid lines are orthogonal.
It is clear that function vfun computed with higher error and we need to improve code for f. On the other hand we can use more simple approach based on PI algorithm.

The result of the reg2.

Comments have been moved to chat; please do not continue the discussion here. — Kuba, Commented Jan 30 at 5:20
(+1) Interesting. I am not really surprised that my FEM model is not really accurate in the corners. But I don't get why vfun is so much off from what it should be. (I will have to think about this at a later time.) — Henrik Schumacher, Commented Jan 30 at 12:55
@HenrikSchumacher I'm also very surprised. I'm guessing there's a routine bug in the code there. — Alex Trounev, Commented Jan 30 at 15:40
@cvgmt Thank you very much. I also calculated these pictures, but I could not yet fully calculate the inverse transformation. — Alex Trounev, Commented Jan 30 at 23:58

cvgmt · Accepted Answer · 2024-01-30 23:30:45Z

Updated

We use schwarz christoffel formula for triangle. $$ \int_{0}^{z} (\zeta-a_1)^{\alpha_1-1} (\zeta-a_1)^{\alpha_2-1}\,\mathrm{d}\zeta$$
To integrate such integral, we set Assumptions -> {Im[z] > 0} since we mapping the upper plane to the triangle.
And we use (z - (u + I*v))/(z - (u - I*v)) mapping the upper plane to the disk and find its inversion.
The final mapping is the disk to triangle.

Clear["Global`*"];
a1 = 0;
a2 = 1;
a3 = ∞;
{a, b, c} = {{2, 1}, {2, 3}, {0, 4}};
(*{a,b,c}=RandomPolygon[3][[1]];*)
{w1, w2, w3} = {a, b, c} . {1, I};
center = Mean[{w1, w2, w3}];
tri = Triangle[{a, b, c}];
{α1, α2, α3} = {VectorAngle[b - a, c - a], VectorAngle[c - b, a - b], 
  VectorAngle[a - c, b - c]}/π;
f[z_] = A + 
   B*Simplify[
     Integrate[(ζ - a1)^(α1 - 1) (ζ - 
          a2)^(α2 - 1), {ζ, 0, z}, 
      Assumptions -> {Im[z] > 0}]];
sol = NSolve[{f[a1] == w1, f[a2] == w2}];
F[z_] = f[z] /. sol[[1]];
{F[a1], F[a2], Limit[F[t], t -> ∞]} == {w1, w2, 
  w3};(*True*)sol1 = 
 z /. FindRoot[F[z] == center, {z, 1}][[1]] // ReIm;
disk2upper[
   w_] = (z /. Solve[(z - (u + I*v))/(z - (u - I*v)) == w, z][[1]]) /. 
   Thread[{u, v} -> sol1];
(*ParametricPlot[ReIm[disk2upper[x+I*y]]//Evaluate,{x,y}∈\
Disk[{0,0},1 - 10^-10],Mesh->30]*)
disk2triangle[x_, y_] = ReIm[F@disk2upper[x + I*y]];
ParametricPlot[
 disk2triangle[x, y], {x, y} ∈ 
  BoundaryDiscretizeRegion@Disk[{0, 0}, 1 - 10^-10], PlotStyle -> None, 
 BoundaryStyle -> None, 
 MeshFunctions -> {Norm@{#3, #4} &, ArcTan @@ {#3, #4} &}, Mesh -> 30,
  Prolog -> {FaceForm[LightBlue], EdgeForm[Cyan], tri}, 
 PlotRange -> RegionBounds[tri], PlotPoints -> 60, MaxRecursion -> 2]

(* the polar coordinate of triangle to disk *)
(* ParametricPlot[{x, y}, {x, y} ∈ Disk[{0, 0}, 1 - 10^-10], 
 MeshFunctions -> {Function[{x, y}, 
    Norm@{Indexed[disk2triangle[x, y], 1], 
      Indexed[disk2triangle[x, y], 2]}], 
   Function[{x, y}, 
    ArcTan @@ {Indexed[disk2triangle[x, y], 1], 
      Indexed[disk2triangle[x, y], 2]}]}, Frame -> False, 
 Axes -> False, Mesh -> {30, 30}]
*)

Original

Not an answer.
Here we draw the polar coordinate orthogrnal lines on Disk[] and map to reg1. We found that the f1 seems not preserve the angle of orthogrnal lines. ( it means 90 Degree)
The f1=f[reg1] where f in the answer of @AlexTrounev https://mathematica.stackexchange.com/a/297116/72111 or the f1 is equal to the f in the answer of @HenrikSchumacher https://mathematica.stackexchange.com/a/178333/72111 when we replace R with reg1.

f1 = f[reg1];
disktoreg1 = 
  ParametricPlot[{x, y}, {x, y} ∈ reg1, 
   MeshFunctions -> {Function[{x, y}, 
      Norm@{Indexed[f1[x, y], 1], Indexed[f1[x, y], 2]}], 
     Function[{x, y}, 
      ArcTan @@ {Indexed[f1[x, y], 1], Indexed[f1[x, y], 2]}]}, 
   Frame -> False, Axes -> False, Mesh -> {30, 30}];
reg1todisk = 
  ParametricPlot[f1[x, y], {x, y} ∈ reg1, 
   MeshFunctions -> {Function[{x, y}, Norm@{x, y}], 
     Function[{x, y}, ArcTan @@ {x, y}]}, Frame -> False, 
   Axes -> False, Mesh -> {30, 30}];
GraphicsRow[{disktoreg1, reg1todisk}]

We check two mesh lines Mesh -> 2.

By this way we can’t see it since you using too rough mesh. — Alex Trounev, Commented Jan 28 at 15:37
@AlexTrounev To check the conformal mapping, we only need to check two orthogonal lines. — cvgmt, Commented Jan 28 at 15:42
@cvgmt: As I see, the angles at the centre of the disk are preserved. — user64494, Commented Jan 28 at 16:26

Stack Exchange Network

How to find the conformal mapping between two regions

2 Answers 2

Updated

Original

Not the answer you're looking for? Browse other questions tagged
graphics
complex
regions
coordinate-transformation
or ask your own question.

Linked

Hot Network Questions

How to find the conformal mapping between two regions

2 Answers 2

Updated

Original

Not the answer you're looking for? Browse other questions tagged graphicscomplexregionscoordinate-transformation or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
graphics
complex
regions
coordinate-transformation
or ask your own question.