Mallas y elementos: FEM

El FEM es probablemente el método numérico más usado en toda la ingeniería — no solo en óptica, sino en mecánica estructural, fluidos, transferencia de calor, acústica. Nació en los años 50 para calcular esfuerzos en alas de avión, y desde los 90 es la herramienta dominante en electromagnetismo en dominio frecuencial. Su fortaleza: mallas que se adaptan a la geometría, refinamiento donde hace falta, y una base matemática sólida que garantiza convergencia.

La idea

En lugar de una rejilla cartesiana rígida (como FDTD), FEM divide el dominio en elementos — triángulos en 2D, tetraedros en 3D — que se adaptan a los contornos del problema. Dentro de cada elemento, la solución se aproxima como una combinación lineal de funciones base simples (típicamente polinomios de bajo orden).

La ecuación de onda vectorial de Helmholtz:

\nabla \times \left(\frac{1}{\mu_r} \nabla \times \mathbf{E}\right) - k_0^2 \varepsilon_r \mathbf{E} = 0

se convierte, vía la formulación débil (multiplicar por una función test e integrar), en un sistema lineal:

\left(\mathbf{K} - k_0^2 \mathbf{M}\right) \mathbf{e} = \mathbf{f}

donde $\mathbf{K}$ es la matriz de rigidez (integrales del rotacional), $\mathbf{M}$ la matriz de masa (integrales de las funciones base ponderadas por $\varepsilon_r$ ), y $\mathbf{e}$ son los coeficientes de la expansión. Ambas matrices son dispersas — cada nodo solo interactúa con sus vecinos directos.

De la formulación fuerte a la débil

La formulación «fuerte» exige que la ecuación de Helmholtz se cumpla punto a punto. La formulación «débil» la relaja: multiplicamos por una función test $\mathbf{w}$ e integramos sobre el dominio $\Omega$ :

\int_\Omega \mathbf{w} \cdot \left[\nabla \times \left(\frac{1}{\mu_r}\nabla \times \mathbf{E}\right) - k_0^2 \varepsilon_r \mathbf{E}\right] d\Omega = 0

Aplicamos integración por partes (el teorema de Green vectorial) al primer término:

\int_\Omega \frac{1}{\mu_r}(\nabla \times \mathbf{w})\cdot(\nabla \times \mathbf{E})\, d\Omega - k_0^2 \int_\Omega \varepsilon_r \, \mathbf{w} \cdot \mathbf{E}\, d\Omega = \oint_{\partial\Omega} \mathbf{w} \cdot \left(\hat{n} \times \frac{1}{\mu_r}\nabla \times \mathbf{E}\right) dS

El término de superficie incorpora las condiciones de contorno (PML, radiación, simetría). El término de volumen, al expandir $\mathbf{E}$ y $\mathbf{w}$ en funciones base $\{\mathbf{N}_i\}$ , da exactamente el sistema $(\mathbf{K} - k_0^2\mathbf{M})\mathbf{e} = \mathbf{f}$ .

La malla

La malla es el corazón de FEM. Una buena malla captura la geometría con precisión, refina donde los campos varían rápido, y mantiene los triángulos bien formados (no muy alargados). Explora cómo se ve una malla para distintas geometrías:

Explorar

Geometría

Densidad 10

Malla de elementos finitos

FEM divide el dominio en triángulos. Cada triángulo lleva una función base lineal. La malla se adapta a la geometría.

Fíjate en la guía de onda con el agujero circular: la malla se adapta al contorno curvo. En FDTD, esa curva se aproximaría con escaleras. En FEM, los triángulos siguen la geometría fielmente. Esta es la ventaja fundamental de las mallas no estructuradas.

Elementos de arista: el truco para el electromagnetismo

En mecánica estructural, FEM usa funciones base nodales — cada nodo almacena un valor escalar. Pero para campos vectoriales como $\mathbf{E}$ , las funciones nodales tienen un problema grave: imponen continuidad de todas las componentes del campo, cuando Maxwell solo requiere continuidad de las componentes tangenciales.

La solución son los elementos de arista (elementos de Nédélec), propuestos en 1980: funciones base asociadas a las aristas del triángulo, no a los nodos. Cada arista almacena la componente tangencial del campo. Esto garantiza automáticamente las condiciones de contorno correctas y elimina las temidas soluciones espurias — modos no físicos que plagaban las implementaciones FEM tempranas.

Ventajas y limitaciones

Mallas adaptativas: puedes poner más elementos donde los campos cambian rápido (cerca de esquinas, interfaces, resonancias) y menos donde son suaves. El refinamiento hp (malla + orden polinómico) da convergencia exponencial.
Matrices dispersas: a diferencia de DDA (matriz densa), FEM produce matrices donde la mayoría de las entradas son cero. Con solucionadores directos (MUMPS, PARDISO), el coste escala como $O(N^{1.5})$ en 3D.
Formulación variacional: la formulación débil garantiza que la solución minimiza un error funcional. Esto da estimaciones de error rigurosas y convergencia demostrable — FEM tiene la base matemática más sólida de todos los métodos numéricos.
Dominios abiertos: FEM necesita truncar el dominio computacional y añadir PML o condiciones de radiación. Esto añade complejidad y grados de libertad adicionales — el siguiente método (BEM) lo resuelve de forma natural.
Dominio frecuencial: cada longitud de onda requiere resolver un sistema lineal independiente. Para respuestas broadband, FDTD es más eficiente (una simulación → todo el espectro).

FEM en la práctica

El software más usado en fotónica es COMSOL Multiphysics (comercial, interfaz gráfica, multifísica). En microondas y antenas: HFSS (Ansys). Específicos para nanofotónica: JCMsuite (hp-FEM, usado extensivamente en metrología y litografía EUV) y FEniCS (open-source, Python, para investigación).

Un cálculo típico: modo guiado de una fibra óptica. La malla tiene ~10.000 triángulos (más finos cerca del núcleo), el sistema resultante es ~30.000 DOF, y se resuelve en segundos. Para un cristal fotónico 3D con defectos, la malla puede tener 10⁶ tetraedros y el sistema 10⁷ DOF — horas en un cluster.

FEM es el método de elección para problemas en dominio frecuencial con geometrías complejas: fibras ópticas, resonadores, cristales fotónicos, guías de onda. La combinación de mallas adaptativas, elementos de arista y convergencia hp lo hace imbatible para problemas estacionarios.

Ejercicios

Ejercicio 1

Observa la malla de la guía de onda con agujero circular en el visualizador de arriba. Cuenta aproximadamente cuántos triángulos hay dentro del agujero vs fuera. ¿Por qué la malla es más fina cerca de la interfaz curva? Si la guía tiene un ancho de 1 μm y cada triángulo tiene un lado de ~50 nm cerca del borde, ¿cuántos elementos por longitud de onda tienes a $\lambda = 1550$ nm?

Solución

La malla es más fina cerca de la interfaz porque el campo cambia rápidamente allí (la discontinuidad de $\varepsilon$ causa un cambio abrupto en la componente normal de $\mathbf{E}$ ). Lejos de la interfaz, el campo varía suavemente y se puede representar con elementos grandes.

A λ = 1550 nm con triángulos de 50 nm: $1550/50 \approx 31$ elementos por λ — más que suficiente (el mínimo práctico para FEM lineal es ~6–10 por λ; con elementos de orden p ≥ 2, bastante menos).

Ejercicio 2

En FEM, cada arista de la malla tiene un grado de libertad (DOF) para elementos de arista de primer orden. Si una malla 2D tiene $N_T$ triángulos, el número de aristas es aproximadamente $1.5 N_T$ (cada triángulo tiene 3 aristas, cada arista se comparte entre 2 triángulos). ¿Cuántos DOF tiene una malla de 10.000 triángulos? ¿Y si usas elementos de segundo orden ( $p = 2$ , dos DOF por arista)?

Solución

Aristas: $\approx 1.5 \times 10\,000 = 15\,000$ . DOF con p=1: 15.000. DOF con p=2: 30.000.

La ventaja de subir el orden (p-refinement) es que la convergencia es $O(h^{p+1})$ , así que duplicar el orden puede ser más eficiente que cuadruplicar el número de triángulos. El coste adicional es relativamente pequeño: las matrices son más grandes pero tienen la misma estructura dispersa.