```python
%pylab inline

from mpl_toolkits.mplot3d import Axes3D
import sympy as sym
```

    Populating the interactive namespace from numpy and matplotlib


Solving 

$$
-\Delta u = f
$$

Example 1 (peak 2d): 

$$
u(x,y) = x(x-1)y(y-1)e^{-100\big( (x-0.5)^2+(y-0.117)^2\big)}
$$

Example 2 (for L-shaped domain, no-rhs, only boundary conditions, in polar coordinates):

$$
u(r,\theta) = r^{2/3} \sin\left(\frac{2\theta+5\pi}{3}\right)
$$

etc.

Play around with your own examples!


```python
from sympy import atan2, sin, cos

def u(x,y):
    return x*(x-1)*y*(y-1)*exp(-100*(x-0.5)**2-100*(y-0.117)**2)

def u1(x,y):
    r = (x*x+y*y)**(1/2)
    theta = atan2(x,y)
    return r**(2/3)*sin(theta*2/3)

def u2(x,y,z):
    return x*(x-1)*y*(y-1)*z*(z-1)*exp(-100*(x-0.5)**2-100*(y-0.117)**2-100*(z-0.331)**2)

def u3(x,y,z):
    return (x*x + y*y + z*z)**(1/4);

def u4(x,y):
    return sin(2*sym.pi*y)*sin(2*sym.pi*x)

def u5(x,y):
    return x*(x-1)*y*(y-1)*(10*exp(-100*(x-0.5)**2-100*(y-0.1)**2)+.1*exp(-100*(x-0.1)**2-100*(y-0.5)**2))

all_2d_us = [u, u1, u4, u5]
all_3d_us = [u2, u3]
```


```python
for a in all_2d_us:
    print(sym.latex(a))
```

    <function u at 0x7fb701623f28>
    <function u1 at 0x7fb701623ea0>
    <function u4 at 0x7fb701623b70>
    <function u5 at 0x7fb701623ae8>



```python
def prm(symu, f):
    print(
        "Exact solution expression     = "+ str(symu).replace('**', '^')+"\n" +
        "Forcing term expression       = "+ str(f).replace('**', '^') + "\n" +
        "Boundary condition expression = "+ str(symu).replace('**', '^')+ "\n" +
        "Problem constants             = "+ "pi:"+str(pi)
    )

def generate_prm(symu):
    f = -symu.diff(x,2)-symu.diff(y,2)
    f = sym.simplify(f)
    prm(symu, f)
```


```python
from sympy import exp
x,y,z = sym.var('x,y,z')

generate_prm(u(x,y))
```

    Exact solution expression     = x*y*(x - 1)*(y - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)
    Forcing term expression       = (x*(x - 1)*(-y*(y - 1)*(200*y - 23.4)^2 + 200*y*(y - 1) + 2*y*(200*y - 23.4) + 2*(y - 1)*(200*y - 23.4) - 2) + y*(y - 1)*(-x*(x - 1)*(200*x - 100.0)^2 + 200*x*(x - 1) + 2*x*(200*x - 100.0) + 2*(x - 1)*(200*x - 100.0) - 2))*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)
    Boundary condition expression = x*y*(x - 1)*(y - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2)
    Problem constants             = pi:3.141592653589793



```python
generate_prm(u1(x,y))
```

    Exact solution expression     = (x^2 + y^2)^0.333333333333333*sin(2*atan2(x, y)/3)
    Forcing term expression       = 0
    Boundary condition expression = (x^2 + y^2)^0.333333333333333*sin(2*atan2(x, y)/3)
    Problem constants             = pi:3.141592653589793



```python
generate_prm(u2(x,y,z))
```

    Exact solution expression     = x*y*z*(x - 1)*(y - 1)*(z - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)
    Forcing term expression       = z*(z - 1)*(x*(x - 1)*(-y*(y - 1)*(200*y - 23.4)^2 + 200*y*(y - 1) + 2*y*(200*y - 23.4) + 2*(y - 1)*(200*y - 23.4) - 2) + y*(y - 1)*(-x*(x - 1)*(200*x - 100.0)^2 + 200*x*(x - 1) + 2*x*(200*x - 100.0) + 2*(x - 1)*(200*x - 100.0) - 2))*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)
    Boundary condition expression = x*y*z*(x - 1)*(y - 1)*(z - 1)*exp(-100*(x - 0.5)^2 - 100*(y - 0.117)^2 - 100*(z - 0.331)^2)
    Problem constants             = pi:3.141592653589793



```python
generate_prm(u3(x,y,z))
```

    Exact solution expression     = (x^2 + y^2 + z^2)^0.25
    Forcing term expression       = (0.75*(x^2 + y^2)*(x^2 + y^2 + z^2)^0.75 - 1.0*(x^2 + y^2 + z^2)^1.75)*(x^2 + y^2 + z^2)^(-2.5)
    Boundary condition expression = (x^2 + y^2 + z^2)^0.25
    Problem constants             = pi:3.141592653589793



```python
generate_prm(u4(x,y))
```

    Exact solution expression     = sin(2*pi*x)*sin(2*pi*y)
    Forcing term expression       = 8*pi^2*sin(2*pi*x)*sin(2*pi*y)
    Boundary condition expression = sin(2*pi*x)*sin(2*pi*y)
    Problem constants             = pi:3.141592653589793



```python
generate_prm(u5(x,y))
```

    Exact solution expression     = x*y*(x - 1)*(y - 1)*(10*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))
    Forcing term expression       = -x*(x - 1)*(y*(y - 1)*(0.1*(200*y - 100.0)^2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + 10*(200*y - 20.0)^2*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 2000*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 20.0*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - 2*y*((20.0*y - 10.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*y - 200.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) - 2*(y - 1)*((20.0*y - 10.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*y - 200.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) + 20*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - y*(y - 1)*(x*(x - 1)*(10*(200*x - 100.0)^2*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*(200*x - 20.0)^2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) - 2000*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) - 20.0*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2)) - 2*x*((20.0*x - 2.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*x - 1000.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) - 2*(x - 1)*((20.0*x - 2.0)*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2) + (2000*x - 1000.0)*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2)) + 20*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.2*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))
    Boundary condition expression = x*y*(x - 1)*(y - 1)*(10*exp(-100*(x - 0.5)^2 - 100*(y - 0.1)^2) + 0.1*exp(-100*(x - 0.1)^2 - 100*(y - 0.5)^2))
    Problem constants             = pi:3.141592653589793



```python

```

c829055c2df0ccef26f69200beeb4cb87e32918a
ipynb
Jupyter Notebook
manufactured_solution.ipynb
dealii-courses/sissa-mhpc-lab-06-iprusak
6a252a81b83277708e7a83b6394a385aebe6f444
MIT
manufactured_solution.ipynb
dealii-courses/sissa-mhpc-lab-06-iprusak
6a252a81b83277708e7a83b6394a385aebe6f444
MIT
manufactured_solution.ipynb
dealii-courses/sissa-mhpc-lab-06-iprusak
6a252a81b83277708e7a83b6394a385aebe6f444
MIT
2021-09-17T08:28:00.000Z
2021-09-17T08:28:00.000Z
Qwen/Qwen-72B
1. YES
2. YES