Ir al contenido principal

Nicolás' blog

Regresando al Método de elementos de frontera

Nicolás Guarín-Zapata

Comentarios

También disponible en:

English

Durante octubre (2017) estuve escribiendo un programa por día para algunos métodos numéricos famosos en Python y Julia. Esto estaba pensado como un ejercicio para aprender Julia. Puedne ver un resumen aquí. Tuve éxito en 30 de los retos, exceptuando por el BEM (Método de elementos de frontera, en inglés), donde no pude encontrar que hice mal ese día. La publicación original está aquí.

Thomas Klimpel encontró el error y me escribió un email describiendo lo que estaba mal. Entonces, estoy creando esta publicación con la implementación correcta del BEM.

El Método de Elementos de Frontera

Queremos rsolver la ecuación

\begin{equation*} \nabla^2 u = -f(x, y)\quad \forall (x, y) \in \Omega\, , \end{equation*}

con

\begin{equation*} u(x) = g(x, y), \quad \forall (x, y)\in \partial \Omega \, . \end{equation*}

Para este método, necesitamos usar la representación integral de la ecuación, que, en este caso, es

\begin{equation*} u(\boldsymbol{\xi}) = \int_{S} [u(\mathbf{x}) F(\mathbf{x}, \boldsymbol{\xi}) - q(\mathbf{x})G(\mathbf{x}, \boldsymbol{\xi})]\mathrm{d}S(\mathbf{x}) + \int_{V} f(\mathbf{x}) G(\mathbf{x}, \boldsymbol{\xi}) \mathrm{d}V(\mathbf{x}) \end{equation*}

con

\begin{equation*} G(\mathbf{x}, \boldsymbol{\xi})= -\frac{1}{2\pi}\ln|\mathbf{x} - \xi| \end{equation*}

y

\begin{equation*} F(\mathbf{x}, \boldsymbol{\xi}) = -\frac{1}{2\pi |\mathbf{x} - \boldsymbol{\xi}|^2} (\mathbf{x} - \boldsymbol{\xi})\cdot\hat{\mathbf{n}} \end{equation*}

Entonces, podemos formar un sistema de ecuaciones

\begin{equation*} [G]\{q\} = [F]\{u\}\, , \end{equation*}

que obtenemos discretizando la frontera. Si tomamos variables constantes en la discretización, las integrales se pueden calcular analíticamente

\begin{equation*} G_{nm} = -\frac{1}{2\pi}\left[r \sin\theta\left(\ln|r| - 1\right) + \theta r\cos\theta\right]^{\theta_B, r_B}_{\theta_A, r_A} \end{equation*}

y

\begin{equation*} F_{nm} = \left[\frac{\theta}{2\pi}\right]^{\theta_B}_{\theta_A} \end{equation*}

para puntos \(n\) y \(m\) en diferentos elementos, y los subíndices \(A, B\) se refieren a los puntos finales del elemento de evaluación. Debemos tener cuidado al evaluar esta expresión ya que \(r_A\) y \(r_B\) pueden ser (casi) cero y hacer que explote. Además, aquí es donde estaba el problema principal ya que olvidé calcular los ángulos para elements que están, en general, no alineados con los ejes vertical u horizontal.

Para los términos en la diagonal las integrales se evalúan como

\begin{equation*} G_{nn} = -\frac{L}{2\pi}\left(\ln\left\vert\frac{L}{2}\right\vert - 1\right) \end{equation*}

y

\begin{equation*} F_{nn} = - \frac{1}{2} \end{equation*}

con \(L\) el tamaño del elemento.

A continuación está el código. Tengan en cuenta que este código funciona para problemas puramente de Dirichlet. Para problemas mixtos Dirichlet-Neumann, las matrices de influencia deben reordenarse para separar las incógnitas de las variables conocidas en lados opuestos de la ecuación.

Pueden descargar los archivos para este proyecto aquí. Se incluye un archivo YML para crear un ambiente de conda con la lista de dependencias. Por ejemplo, se usa la versión 3.0 de Meshio.

import numpy as np
from numpy import log, arctan2, pi, mean, arctan
from numpy.linalg import norm, solve
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import meshio


  def assem(coords, elems):
      """Assembly matrices for the BEM problem

      Parameters
      ----------
      coords : ndarray, float
          Coordinates for the nodes.
      elems : ndarray, int
          Connectivity for the elements.

      Returns
      -------
      Gmat : ndarray, float
          Influence matrix for the flow.
      Fmat : ndarray, float
          Influence matrix for primary variable.
      """
      nelems = elems.shape[0]
      Gmat = np.zeros((nelems, nelems))
      Fmat = np.zeros((nelems, nelems))
      for ev_cont, elem1 in enumerate(elems):
          for col_cont, elem2 in enumerate(elems):
              pt_col = mean(coords[elem2], axis=0)
              if ev_cont == col_cont:
                  L = norm(coords[elem1[1]] - coords[elem1[0]])
                  Gmat[ev_cont, ev_cont] = - L/(2*pi)*(log(L/2) - 1)
                  Fmat[ev_cont, ev_cont] = - 0.5
              else:
                  Gij, Fij = influence_coeff(elem1, coords, pt_col)
                  Gmat[ev_cont, col_cont] = Gij
                  Fmat[ev_cont, col_cont] = Fij
      return Gmat, Fmat


  def influence_coeff(elem, coords, pt_col):
      """Compute influence coefficients

      Parameters
      ----------
      elems : ndarray, int
          Connectivity for the elements.
      coords : ndarray, float
          Coordinates for the nodes.
      pt_col : ndarray
          Coordinates of the colocation point.

      Returns
      -------
      G_coeff : float
          Influence coefficient for flows.
      F_coeff : float
          Influence coefficient for primary variable.
      """
      dcos = coords[elem[1]] - coords[elem[0]]
      dcos = dcos / norm(dcos)
      rotmat = np.array([[dcos[1], -dcos[0]],
                      [dcos[0], dcos[1]]])
      r_A = rotmat.dot(coords[elem[0]] - pt_col)
      r_B = rotmat.dot(coords[elem[1]] - pt_col)
      theta_A = arctan2(r_A[1], r_A[0])
      theta_B = arctan2(r_B[1], r_B[0])
      if norm(r_A) <= 1e-6:
          G_coeff = r_B[1]*(log(norm(r_B)) - 1) + theta_B*r_B[0]
      elif norm(r_B) <= 1e-6:
          G_coeff = -(r_A[1]*(log(norm(r_A)) - 1) + theta_A*r_A[0])
      else:
          G_coeff = r_B[1]*(log(norm(r_B)) - 1) + theta_B*r_B[0] -\
                  (r_A[1]*(log(norm(r_A)) - 1) + theta_A*r_A[0])
      F_coeff = theta_B - theta_A
      return -G_coeff/(2*pi), F_coeff/(2*pi)


  def eval_sol(ev_coords, coords, elems, u_boundary, q_boundary):
      """Evaluate the solution in a set of points

      Parameters
      ----------
      ev_coords : ndarray, float
          Coordinates of the evaluation points.
      coords : ndarray, float
          Coordinates for the nodes.
      elems : ndarray, int
          Connectivity for the elements.
      u_boundary : ndarray, float
          Primary variable in the nodes.
      q_boundary : [type]
          Flows in the nodes.

      Returns
      -------
      solution : ndarray, float
          Solution evaluated in the given points.
      """
      npts = ev_coords.shape[0]
      solution = np.zeros(npts)
      for k in range(npts):
          for ev_cont, elem in enumerate(elems):
              pt_col = ev_coords[k]
              G, F = influence_coeff(elem, coords, pt_col)
              solution[k] += u_boundary[ev_cont]*F - q_boundary[ev_cont]*G
      return solution


#%% Simulation
mesh = meshio.read("disk.msh")
elems = mesh.cells["line"]
bound_nodes = list(set(elems.flatten()))
coords = mesh.points[bound_nodes, :2]
x, y = coords.T
x_m, y_m = 0.5*(coords[elems[:, 0]] + coords[elems[:, 1]]).T
theta = np.arctan2(y_m, x_m)
u_boundary = 3*np.cos(6*theta)


#%% Assembly
Gmat, Fmat = assem(coords, elems)

#%% Solution
q_boundary = solve(Gmat, Fmat.dot(u_boundary))

#%% Evaluation
ev_coords =  mesh.points[:, :2]
ev_x, ev_y = ev_coords.T
solution = eval_sol(ev_coords, coords, elems, u_boundary, q_boundary)

#%% Visualization
tris = mesh.cells["triangle"]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_trisurf(ev_x, ev_y, solution, cmap="RdYlBu", lw=0.3,
                edgecolor="#3c3c3c")
plt.xticks([])
plt.yticks([])
ax.set_zticks([])
plt.savefig("bem_solution.png", bbox_inches="tight", transparent=True,
            dpi=300)

En este caso obtenemos el siguiente resultado.

Solución de la ecuación diferencial usando BEM.

Comentarios

Comments powered by Disqus