# Ecuación de Poisson

INTRODUCCION

## La ecuación diferencial parcial más sencilla

Solucionaremos la ecuación de Poisson, el problema "Hola Mundo" de la computación científica:  
    Encontrar la incógnita $u$, conociendo el término fuente $f$, tal que:
    \begin{align}
      -\triangle u &= f, \ && \text{en} \: \Omega,\\
      u &= u_0, && \text{sobre} \: \partial\Omega.
    \end{align}


La ecuación modela distintos fenómenos difusivos:

- Conducción de calor, electrostática, flujos no viscosos, ondas.
- Como parte de métodos de solución de EDPs para configuraciones no-clásicas.

## La receta numérica para aproximar su solución

Aproximar numéricamente el problema por métodos numéricos variacionales consta de los siguiente pasos:
- Plantear la formulación variacional del problema.
- Plantear el espacio funcional de aproximación numérica: Galerkin Continuo, Galerkin Discontinuo, Colocación, etc.
- Ensamblar el sistema discreto de ecuaciones, incluyendo las condiciones de contorno.
FEniCS nos reduce este proceso al primer paso!

## Formulación variacional

La formulación variacional la definimos como la integral sobre el dominio $\Omega$ del producto interno entre la ecuación de Poisson y una función vectorial de test $\boldsymbol{v}\in \hat{V}$, donde $\hat{V}$ es el espacio de la funcion de test.

\begin{align*} 
-\int_\Omega (\triangle u) v \text{d} x &= \int_\Omega fv \text{d} x, && \text{en} \: \Omega.
\end{align*}

Como $\triangle u$ contiene derivadas de segundo orden sobre $u$, debilitamos la ecuación integrando por partes:

\begin{align*}
-\int_\Omega (\triangle u) v \text{d} x &= \int_\Omega \nabla u \cdot \nabla v \text{d} x - \int_{\partial \Omega}\frac{\partial u}{\partial n} v \text{d} s, && \text{en} \: \Omega.
\end{align*}
donde $\boldsymbol{n}$ es vector normal en el contorno.

$\frac{\partial u}{\partial n}$ es conocido como el flujo en el contorno $\partial \Omega$ y puede ser usado para prescribir una condicion de contorno de tipo Neumann o Robin de la forma $\frac{\partial u}{\partial n} = t$.

Suponiendo que no hay flujos en los contornos del dominio, la formulación variacional está dada por
    \begin{align*}
         \int_\Omega \nabla u \cdot \nabla v \text{d} x&= \int_\Omega fv \text{d} x, && \text{en} \: \Omega.\\
    \end{align*}
    
### Definición de los espacios discretos de interpolación

Definición de un espacio discreto de aproximación contenido dentro del espacio continuo. Uno para el espacio de la función de test:
    \begin{align*}
        \widetilde{V}_h\subset \widetilde{V}:=\left\{v\in H(\Omega): v=0 \:\text{sobre} \: \partial\Omega\right\}, 
    \end{align*}
    y otra para la función de la incógnita:
    \begin{align*}
        V_h\subset V:=\left\{v\in H(\Omega): v=u_0 \:\text{sobre} \: \partial\Omega\right\}.
    \end{align*}
Tal que, el problema discreto es encontrar $u_h\in V_h\subset V$, tal que:
    \begin{align*}
     \int_\Omega \nabla u_h \cdot \nabla v_h \text{d} x&= \int_\Omega fv_h \text{d} x, && \text{en} \: \Omega.\\
    \end{align*}
    para todo $v_h\in \widetilde{V}_h\subset \widetilde{V}$.


### Notación variacional en FEniCS
Este tipo de ecuaciones variacionales discretas se puede escribir usando la siguiente notación: Encontrar $u_h\in V_h\subset V$, tal que:
    \begin{align*}
     a(u_h,v_h)&=L(v_h),
    \end{align*}
 para todo $v_h\in \widetilde{V}_h\subset \widetilde{V}$.
 
 Donde $a(u,v)$ es una forma variacional bilineal (lineal en cada argumento) y $L(v)$ es una forma lineal:
    \begin{align*}
     a(u_h,v_h)&=\int_\Omega \nabla u_h \cdot \nabla v_h \text{d} x,\\
     L(v_h)&= \int_\Omega fv_h \text{d} x.\\
     \end{align*}
     
  ## Sistema discreto de ecuaciones lineales

Introduciendo las funciones de test y de interpolación de la incógnita como polinomios a trozos:
\begin{align*}
        v_h=&\sum_i^n \phi_i \subset \widetilde{V}_h\subset \widetilde{V}, &&
        u_h=&\sum_j^n U_j\phi_j \subset V_h\subset V.
    \end{align*}
La ecuación discreta resulta en:
\begin{align*}
    \int_\Omega \partial_k \sum_j^n U_j\phi_j \partial_k \sum_i^n \phi_i \text{d} x =& \int_\Omega f\sum_i \phi_i \text{d} x.
\end{align*}
O mejor:
\begin{align*}
    \sum_i^n\sum_j^n\int_{\Omega^e} U_j \partial_k\phi_j \partial_k \phi_i \text{d} x =& \sum_i^n \int_{\Omega^e} f \phi_i \text{d} x.
\end{align*}


Deshaciendo las sumatorias y el álgebra lineal:
\begin{align*}
        \begin{bmatrix}
a_{11} & a_{12} & a_{13} & \dots & a_{1n} \\
a_{21} & a_{22} & a_{23} & \dots & a_{2n} \\
\dots  & \dots  & \dots  & a_{ij} & \dots  \\
a_{n1} & a_{n2} & a_{n3} & \dots & a_{nn} 
\end{bmatrix}
\begin{bmatrix}
U_1 \\ U_2 \\ \dots \\ U_n 
\end{bmatrix}
=&
\begin{bmatrix}
L_{1} \\ L_{2} \\ \dots \\ L_{3}
\end{bmatrix},
\end{align*}
donde $U_j$ son las incógnitas nodales, $a_{ij}$ viene dada por una integración local sobre cada elemento $a_{ij} = \int_{\Omega^e} U_j \partial_k\phi_j \partial_k \phi_i \text{d} x$ y $L_i$ está dada por $L_i = \int_{\Omega^e} f \phi_i \text{d} x$.

Esto es lo mismo que un sistema lineal del tipo $\mathbb{A}\mathbb{U}=\mathbb{L}$, cuya resolución está dada por:

\begin{align*}
    \mathbb{U}=&\mathbb{A}^{-1}\mathbb{L}
\end{align*}

## Implementación en FEniCS

Modelamos un dominio bidimensional cuadrador con condiciones de Dirichlet en sus contornos y un término fuente constante.




### El programa completo


```python
from dolfin import *
import matplotlib.pyplot as plt

# Malla
mesh = UnitSquareMesh(8,8)
# Espacio de funciones para la incognita y test
V= FunctionSpace(mesh,"Lagrange",1)
# Expresion para las condiciones de contorno
u0=Expression("1+ x[0]*x[1] + 2*x[1]*x[1]",degree=2)
# Definicion de las condiciones de contorno
bc = DirichletBC ( V , u0 , "on_boundary")
# Fuente
f=Constant(-6.0)

# Funcion de trial
u= TrialFunction(V)
# Funcion de test
v= TestFunction(V)
# Forma bilineal del lado izquierdo
a = inner(grad(u), grad(v))*dx
# Forma lineal del lado derecho
L=f*v*dx
# Funcion donde se almacena el resultado
u = Function(V)
# Solucion del sistema algebraico
solve(a==L, u, bc)

#Impresion usando paraview
vtkfile = File('poisson_solution.pvd')
vtkfile << u
```

    Calling FFC just-in-time (JIT) compiler, this may take some time.
    Calling FFC just-in-time (JIT) compiler, this may take some time.


### Importar paquetes

Importamos los paquetes de solucion (dolfin), de mallado (mshr), numpy para realizar operaciones numericas y matplotlib para ver los resultados


```python
from dolfin import *
import matplotlib.pyplot as plt
```

### Dominio geometrico y mallado

Definimos el dominio $\Omega$ y generamos la malla directamente usando FEniCS:


```python
mesh = UnitSquareMesh(8,8)
```

Las posibilidades que existen de generación de dominios con la librería **mshr** son inmensas.

###  Definición de los espacios funcionales

FEniCS tienen programada internamente una gran cantidad de familias de elementos finitos, incluyendo: Lagrangeanos, Curl, Div, Discontinuos, y la posibilidad de realizar distintas combinaciones de espacios funcionales entre ellos.



Una vez creada la malla, podemos crear un espacio funcional de elementos finitos V sobre dicho dominio geométrico discreto:


```python
V= FunctionSpace(mesh,"Lagrange",1)
```

El segundo argumento especifica el tipo de elemento. El tipo de elemento aquí es "Lagrange", lo que implica la familia de elementos estándar de Lagrange. También puede utilizar 'P' para especificar este tipo de elemento. FEniCS soporta todas las familias de elementos simplex y la notación definida en la Tabla Periódica de Elementos Finitos.

El tercer argumento 1 especifica el grado del elemento finito. En este caso, el estándar 𝖯1 elemento lineal de Lagrange, que es un triángulo con nodos en los tres vértices. Algunos se refieren a este elemento como el "triángulo lineal". La solución calculada $u$ será continua en los elementos y variará linealmente en $(x,y)$ dentro de cada elemento. Las aproximaciones polinomiales de mayor grado sobre cada celda se obtienen trivialmente aumentando el tercer parámetro. Cambiar el segundo parámetro a 'DP' crea un espacio funcional para métodos de Galerkin discontinuos.

### Definición de las funciones de trial y test

En matemáticas, distinguimos entre los espacios de interpolación de la incógnita $V$ y el de prueba o test $\tilde{V}$. La única diferencia en el problema actual son las condiciones de contorno. En FEniCS no especificamos las condiciones de contorno como parte del espacio funcional, por lo que es suficiente trabajar con un espacio común $V$ para las distintas funciones:


```python
u = TrialFunction(V)
v = TestFunction(V)
```

### Definición de condiciones de contorno tipo Dirichlet

El siguiente paso es especificar la condición de contorno $u= u_D$ sobre $\partial\Omega$. Esto se hace muy fácilmente con la siguiente rutina:


```python
bc = DirichletBC(V, u0, "on_boundary")
```

donde $u_D$ es una expresión que define los valores de la solución en el contorno del dominio computacional. Dicho valor se puede determinar por una función (u objeto) definida sobre todo el dominio espacial, pero evaluada únicamente en el contorno.

La clase **DirichletBC** sirve para imponer condiciones de contorno en cada espacio funcional (incógnita).

Define que la función de la incógnita $u$ definida en el espacio $V$ debe igualarse a la expresión $u_D$ evaluada en la región "on_boundary".

Las condiciones de contorno del tipo $u = uD$ se conocen como condiciones de Dirichlet. Para el presente método de elementos finitos para el problema de Poisson, también se denominan condiciones de contorno esenciales, ya que deben imponerse explícitamente como parte del espacio funcional de la incógnita (en contraste con las que pueden definirse implícitamente como parte de la formulación variacional). Naturalmente, la clase FEniCS utilizada para definir las condiciones de contorno de Dirichlet se denomina DirichletBC. Esa clase instancia un objeto bc que se usa más adelante para ensamblar el sistema de ecuaciones lineales.

Ahora bien, la variable $u_D$ se refiere a un objeto **Expression**, que se utiliza para representar una función matemática. La construcción típica de ese tipo de objeto se da con la siguiente rutina:


```python
u0 = Expression(formula, degree=1)
```

donde **formula** es una cadena que contiene una expresión matemática. La fórmula debe escribirse con sintaxis C ++ y se convierte automáticamente en una función C ++ compilada y eficiente.

La expresión puede depender de las variables x[0] y x[1] correspondientes a las coordenadas $(x,y)$. En 3D, la expresión también puede depender de la variable x[2] correspondiente a la coordenada $z$. Este problema introductorio supone una expresión de la forma $u_D(x,y) = 1 + x2 + 2y^2$. Entonces, la expresión se puede escribir como 1 + x [0] * x [0] + 2 * x [1] * x [1]:


```python
u0=Expression("1+ x[0]*x[1] + 2*x[1]*x[1]", degree=2)
```

Establecemos el grado **degree** de la fórmula como 2 para que $u_D$ pueda representar la solución cuadrática exacta de nuestro problema de prueba.

La clase **Expression** es muy flexible y sirve para crear expresiones complicadas definidas por el usuario. Funciones espaciales, temporales u otras.

Por ejemplo, la fuente $f=-6$ puede ser definida como:


```python
f = Expression("-6",degree=0)
```

O más eficientemente con la clase Función Constante (en espacio y tiempo):


```python
f = Constant(-6.0)
```

Lo único que nos falta es definir la región sobre la cual se evalúa la expresión para $u_D$. Esto se realiza con el argumento "on_boundary". Otras opciones para determinar fácilmente la región de imposición de las condiciones de Dirichlet son:

- "on\_boundary"    Todo el contorno del dominio
- "near(x[0],0.0)"      En la región del contorno x=0.
- "near(x[0],0.0) || near(x[1],1.0)"    En la región del contorno x=0 o y=1.

### Definición del problema variacional

Ahora tenemos todos los ingredientes que necesitamos para definir el problema variacional:


```python
a = inner(grad(u), grad(v))*dx
L = f*v*dx
```

En esencia, estas dos líneas especifican el problema EDP a resolver. Se entiende la estrecha correspondencia entre la sintaxis de Python y las fórmulas matemáticas $\nabla u \cdot \nabla v dx$ y $fvdx$. Esta es una fortaleza clave de FEniCS: las fórmulas en la formulación variacional se traducen directamente a un código Python muy similar, una característica que facilita la formulación y resolución de problemas EDP complicados.

### Ensamblaje y solución del problema variacional

FEniCS realiza el ensamblaje del sistema lineal (impone las condiciones de contorno) y llama a uno de los solucionadores externos (según el problema).


```python
# Compute solution
u = Function(V)
# Solucion del sistema algebraico
solve(a==L, u, bc)
```

Se usa la función $u$ tanto como función trial como para guardar la solución que entrega el solucionador. 

## Post-procesamiento

### Visualización de la solución

Posproceso, visualización y salida de datos puede realizarse:
- Interno: Numpy, Matplotlib, etc.
- Externo: .pvd, .xml, .xdmf (reestablecer casos).

Por ejemplo, para escribir un archivo que pueda ser visualizado en Paraview, se escriben las dos líneas siguientes:


```python
vtkfile = File("poisson_solution.pvd")
vtkfile << u
```

Una para crear el archivo que guarda los datos.
La otra línea para salvar la información dentro del archivo creado.

e9c569b3714f7d7f191337a1d09c08f9f1597756
ipynb
Jupyter Notebook
poisson.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
poisson.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
poisson.ipynb
ciaid-colombia/Taller-FEniCS
629b26c4a756d1017e6b409e5ed700498a9edacd
MIT
2021-04-24T14:08:38.000Z
2021-04-24T14:08:38.000Z
Qwen/Qwen-72B
1. YES
2. YES