# How to find the conformal mapping between two regions

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.

• Look in mathoverflow.net/questions/314189/… for info. Commented Apr 4, 2023 at 5:42
• Isn't this answer mathematica.stackexchange.com/a/178333/9469 completely solves your problem? 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. 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

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["NDSolveFEM"];
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["NDSolveFEM"];
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 = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y}}];
sd = NDSolveSolutionData[{"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"];

(*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 =
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["NDSolveFEM"];
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.) Commented Jan 30 at 12:55
• @HenrikSchumacher I'm also very surprised. I'm guessing there's a routine bug in the code there. Commented Jan 30 at 15:40
• Add the picture of the reg2. Commented Jan 30 at 23:35
• @cvgmt Thank you very much. I also calculated these pictures, but I could not yet fully calculate the inverse transformation. Commented Jan 30 at 23:58

## 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`.

• Please indicate what the right angle you mean. Commented Jan 28 at 15:35
• By this way we can’t see it since you using too rough mesh. Commented Jan 28 at 15:37
• @AlexTrounev To check the conformal mapping, we only need to check two orthogonal lines. Commented Jan 28 at 15:42
• @cvgmt We need to check two orthogonal lines in a point ;) Commented Jan 28 at 16:07
• @cvgmt: As I see, the angles at the centre of the disk are preserved. Commented Jan 28 at 16:26