2
$\begingroup$

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.

enter image description here

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

enter image description here

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

enter image description here

and the coefficients

enter image description here

The left and right boundary conditions discretized

enter image description here

enter image description here

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 'My results'

And the expected results

enter image description here

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 
            
    ```
$\endgroup$
1
  • 1
    $\begingroup$ Why do you use the analytical expression for the transformation? This complicates the whole stuff. Simply approximate the metric terms analog to $u$. $\endgroup$
    – ConvexHull
    Dec 9, 2021 at 16:58

2 Answers 2

1
$\begingroup$

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

\begin{equation} \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}. \end{equation}

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

$\endgroup$
0
$\begingroup$

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

\begin{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} \end{equation}

After factoring this is the same \begin{equation} \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} \end{equation}

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}

$\endgroup$
1
  • $\begingroup$ 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}$$ $\endgroup$
    – bc_eng
    Dec 9, 2021 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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