# Possible to use Iterative FD methods to solve a transformed non square domain [matlab]?

For the 2-D Poisson equation $$-(u_{xx}+u_{yy}) = f \ \ \text{where} f = 1$$

For boundary conditions $$\frac{\partial u}{\partial n} = 0 \ \text{on AB and AD}$$ $$u = 0 \ \ \ \text{on BC and CD no-slip condition}$$

I know how to solve this using both a Gauss Siedel and Jacobian method for normal problems however, what if the domain is an irregular shape such as.

I've solved this in class using a transformation of variables and using the Matrix method but it was a pain in the ass

The transformation is $$x = m\xi \\ y = h \eta+s \xi - s\eta \xi$$ Where $$h = 1 s = 0.5 m = 0.866$$ I can post the full math for transforming if interested but the transformed poisson equation is

$$-1/J^2 (au_{\xi \xi} - 2bu_{\xi\eta} + cu_{\eta \eta} + du_{\eta}+eu_\xi)$$ $$\text{Where} \ \ u_\xi \text{ is discretized like } u_x$$ Full Discretization scheme is like

and the coefficients

The left and right boundary conditions discretized

I put together a jacobi scheme like this

$$u^{k+1}_{ij} = \frac{1}{2*h_{\xi\eta}} *\{a\Delta \eta^2(u_{i+1,j}+u_{i-1,j})+c\Delta \xi^2 (u_{i,j+1}+u_{i,j-1}) + \\ [2d\Delta \xi(u_{i,j+1}-u_{i,j-1}) - 2b(u_{i+1,j+1}+u_{i-1,j-1}-u_{i-1,j+1}-u_{u+1,j-1})]\frac{\Delta \xi \Delta \eta}{4} \\ + J^2f_{i,j} \Delta xi^2 \Delta \eta^2 \}$$ $$h_{ \xi \eta} = \Delta \xi^2 + \Delta \eta^2$$

I've iterated about 1000 times but my results don't align with the expected

And the expected results

I think I maybe having issues implementing the boundary conditions so here's my code, but if it's not possible to do this iteratively then just let me know. THANK YOU

%% Initiate Variables
N = 21;

l = 3;
s = 0.5;
h = 1;

m = sqrt((1/2*l-h+s)^2-s^2);
x = linspace(0,m,N);
y = linspace(0,h,N);

dx = x(2)-x(1); dy = y(2)-y(1);

%% Map coordinates to the xi and nu

xi = x/m;
nu = (m.*y-s.*x)./(m*h-s.*x);

dxi = xi(2)-xi(1); %these should be the same i think
dnu = nu(2)-nu(1);

hx = dxi.^2;
hn = dnu.^2;
hh = dxi*dnu; %multipy non squared
hhxn = hn*hx;
hxn = hx+hn; %adding the squres
w = 2/3;
%% Solve the poisson equation I guess ?

u = zeros(N,N);
uj = u;
f = u;

J = m*(h-s*xi);
a = (h-s.*xi).^2;
b = s*(1-nu).*(h-s.*xi);
c = m^2+s^2*(1-nu).^2;
e = 0;
af = 0;
beta = 2*s^2*(1-nu).*(h-s.*xi);
d = -2*s.^2*(1-nu);

gam = (s*(1-nu))./(h-s*xi);

f(1:N,1:N) = 1;

%% Jacobi for transformed index

for k = 1:100

uo = u;

for j = 2:N-1
for i= 2:N-1

dudxi2 = uo(i+1,j) + uo(i-1,j);
dudn2 = uo(i,j+1)+uo(i,j-1);

dudx = uo(i+1,j) - uo(i-1,j);
dudn = uo(i,j+1)-uo(i,j-1);

dudxn = uo(i+1,j+1)+uo(i-1,j-1)+uo(i-1,j+1)-uo(i+1,j-1);

u(i,j) = ( 1/(2*hxn)*(a(i)*hn*dudxi2 + c(j)*hx*dudn2 + ...
(2*d(j)*dxi*dudn  + 2*b(i)*dudxn)*(hh/4)) ...
+ (J(i)^2*f(i,j)*hhxn)/(hxn));

end

end

end

%% Boundary conditions

% for k = 1:100
%     uo = u;
% for j = 2:N-1
%
%     i = 1;
%     k1 = dxi/dnu*gam(j)*(uo(i,j+1) - uo(i,j-1))-uo(i+1,j);
%
%     dudxi2 = uo(i+1,j) + uo(i-1,j);
%             dudn2 = uo(i,j-1)+uo(i,j-1);
%
%             dudx = uo(i+1,j) - k1;
%             dudn = 0;
%
%             dudxn = uo(i+1,j-1)+uo(i-1,j-1)+uo(i-1,j-1)-uo(i+1,j-1);
%
%             u(i,j) = 1/(2*hxn)*(a(i)*hn*dudxi2 + c(j)*hx*dudn2 + ...
%                 (2*d(j)*dxi*dudn  + 2*b(i)*dudxn)*(hh/4)) ...
%                 + (J(i,j)^2*f(i,j)*hhxn)/(hxn);
%
% end
% end
%
% for k = 1:100
%     uo = u;
%     for i = 2:N-1
%         j = N;
%         dudxi2 = uo(i+1,j) + uo(i-1,j);
%             dudn2 = uo(i,j-1)+uo(i,j-1);
%
%             dudx = uo(i+1,j) - uo(i-1,j);
%             dudn = 0;
%
%             dudxn = uo(i+1,j-1)+uo(i-1,j-1)+uo(i-1,j-1)-uo(i+1,j-1);
%
%             u(i,j) = 1/(2*hxn)*(a(i)*hn*dudxi2 + c(j)*hx*dudn2 + ...
%                 (2*d(j)*dxi*dudn  + 2*b(i)*dudxn)*(hh/4)) ...
%                 + (J(i,j)^2*f(i,j)*hhxn)/(hxn);
%     end
% end
%

%% plot

figure()
contour(x,y,u)

for i = 1:N
for j = 1:N
xi2 = (i-1)*dx;
nu2 = (j-1)*dy;
X(i,j) = m*xi2;
Y(i,j) = h*nu2 + s*xi2 - s*nu2*xi2;
end
end
umax = max(max(max(abs(u))))

figure()
contour(X,Y,u');

figure,
mesh(X,Y,0*X,0*Y);
view([0,0,1]);
axis equal

figure,
patch([0,m,m,0],[0,s,h,h],-ones(1,4),0,'facecolor',[0.8,.8,.8]);
hold on
[ccc,fff] = contour(X,Y,abs(u),(0.02:0.02:max(max(max(abs(u)))))');
clabel(ccc,fff)
axis equal
hold off

$$$$
`
• Why do you use the analytical expression for the transformation? This complicates the whole stuff. Simply approximate the metric terms analog to $u$. Dec 9, 2021 at 16:58

A hint for you:

Starting from the Poisson equation

\begin{align} u_{xx} + u_{yy} &= f \quad \text{in} ~~ \Omega , \\ u & = 0 \quad \text{in} ~~ \partial \Omega ,\\ \frac{\partial u}{\partial n} &= 0 \quad \text{in} ~~ \partial \Omega. \end{align}

You can use following general relation

$$$$\underbrace{\begin{pmatrix} x_{\xi} & y_{\xi} & 0 & 0 & 0 \\ x_{\eta} & y_{\eta} & 0 & 0 & 0 \\ x_{\xi\xi} & y_{\xi\xi} & x_{\xi}^2 & y_{\xi}^2 & 2 x_{\xi} y_{\xi} \\ x_{\eta\eta} & y_{\eta\eta} & x_{\eta}^2 & y_{\eta}^2 & 2 x_{\eta} y_{\eta} \\ x_{\xi\eta} & y_{\xi\eta} & x_{\eta} x_{\xi} & y_{\eta} y_{\xi} & x_{\eta} y_{\xi} + x_{\xi} y_{\eta} \end{pmatrix}}_{\underline{\underline{M}}} \cdot \begin{pmatrix} \partial_{x} \\ \partial_{y} \\ \partial_{xx}^2 \\ \partial_{yy}^2 \\ \partial_{x}\partial_{y} \\ \end{pmatrix} = \begin{pmatrix} \partial_{\xi} \\ \partial_{\eta} \\ \partial_{\xi\xi}^2 \\ \partial_{\eta\eta}^2 \\ \partial_{\xi}\partial_{\eta} \end{pmatrix}.$$$$

Do not use the analytical expression for the transformation. Instead discretisize the Metric terms

\begin{align} x_{\xi~(i,j)} \approx \frac{x_{i+1,j} - x_{i-1,j}}{2\Delta\xi}, \quad\quad & y_{\xi~(i,j)} \approx \frac{y_{i+1,j} - y_{i-1,j}}{2\Delta\xi} \\ x_{\eta~(i,j)} \approx \frac{x_{i,j+1} - x_{i,j-1}}{2\Delta\eta}, \quad\quad & y_{\eta~(i,j)} \approx \frac{y_{i,j+1} - y_{i,j-1}}{2\Delta\eta} \\ \dots \end{align}

Use one-side stencils at the boundaries.

Now simply invert the matrix with $$\underline{\underline{M}}^{-1} = \underline{\underline{W}}$$

\begin{align} \underline{\underline{M}} \cdot \underline{\partial_{\boldsymbol{x}}} &= \underline{\partial_{\boldsymbol{\xi}}} ,\\ \underline{\partial_{\boldsymbol{x}}} &= \underline{\underline{W}} \cdot \underline{\partial_{\boldsymbol{\xi}}}, \end{align}

resulting in

$$\begin{eqnarray} a u_{\xi\xi} + b u_{\xi\eta} + c u_{\eta\eta} + d u_{\xi} + e u_{\eta} = f, \label{eq:dgl_trafo} \end{eqnarray}$$

where $$a$$, $$b$$, $$c$$, $$d$$, $$e$$ are point wise values which can be used directly on each DOF $$(i,j)$$. To build up $$a$$, $$b$$, $$c$$, $$d$$, $$e$$ you only need column 3 and 4 of $$\underline{\underline{W}}$$, guess why?

Works on arbitrary meshes!

The rest is quite simple and similar to the Cartesian case. Simply discretisize $$u$$ and use a Gauss-Seidel or SOR method.

Regards

To directly answer your question. Yes, it is possible, but it gets complicated as one needs to keep track of coordinate changes, and the local metric. I have not checked in detail if the analytical formulation of the Laplacian you have is correct for your coordinate system. I did check the boundary condition, and it seems to be correct.

I notice in your code you compute your grid in with constant $$dx$$ and $$dy$$ and then, it seems, assume $$d\xi$$ and $$d\eta$$ are also constant. In one or the other coordinate systems they are not having constant spacing.

I did calculate the boundary conditions and I get the following along $$AC$$ which I think is the same as your discretized version. (Substitute $$\xi=x/m$$ to compare with your equation).

$$$$\frac{\partial u}{\partial n} = \nabla u\cdot (-\mathbf{e}_{x}) = - \left(- \frac{s}{m \left(h - \frac{s x}{m}\right)} + \frac{s \left(y - \frac{s x}{m}\right)}{m \left(h - \frac{s x}{m}\right)^{2}}\right) \frac{\partial}{\partial \eta} u{\left(\xi,\eta \right)} - \frac{\frac{\partial}{\partial \xi} u{\left(\xi,\eta \right)}}{m}$$$$

After factoring this is the same $$$$\frac{\partial u}{\partial n} = \frac{s(1-\eta)}{m(h-s\xi)}\frac{\partial}{\partial \eta} u{\left(\xi,\eta \right)} - \frac{\frac{\partial}{\partial \xi} u{\left(\xi,\eta \right)}}{m}$$$$

The following is not an essential problem. I think your coordinate change does not work the way you think it does. I'm not sure exactly the derivation in your mind, but at I think $$\eta=0$$ the idea is that $$\xi$$ parameterizes the lower boundary $$BC$$. Then similarly as you move up the figure, $$\xi$$ parameterizes straight curves from $$AB$$ to $$CD$$, but these are not all of the same slope if "evenly" distributed, clearly for $$\eta=1$$ the curve $$AC$$ is parallel to the $$x$$-axis while for $$\eta=0$$ $$BC$$ is not.

For the coordinate change you give, $$\frac{\partial x}{\partial \xi} = m$$ for all values of $$\eta$$, but if the intention is as above this should also depend on $$\eta$$. If you think of the slope of these lines as having a smoothly varying angle $$\theta(\eta)$$ such that $$\tan \theta(0) = \frac{s}{m}$$ and $$\tan \theta(1) = 0$$ then one should get the coordinate transform

\begin{aligned} x &= \xi \cos \theta(\eta) \\ y &= h\eta + m\xi\sin \theta(\eta) \end{aligned}

The trignometric value may be found by assuming the base remains $$m$$ for all $$\eta$$ while the height gained across the domain decreases linearly with $$\eta$$ as $$s(1-\eta)$$ gives a hypotenuse of $$h^2=(s(1-\eta))^2 + m^2$$ then $$\cos \theta = \frac{m}{h}$$ and $$\sin \theta = \frac{s(1-\eta)}{h}$$.

Then the transformation would be: \begin{aligned} x&= \frac{m \xi}{\left(m^{2} + s^{2} \left(1 - \eta\right)^{2}\right)^{0.5}} \\ y&=\eta h + \frac{m s \xi \left(1 - \eta\right)}{\left(m^{2} + s^{2} \left(1 - \eta\right)^{2}\right)^{0.5}} \end{aligned}

• No my coordinate change was taken straight out of a textbook $$x \text{ is only a function of } \xi \text{ where as } y \text{ is a function of both}$$ Dec 9, 2021 at 16:19