Archivo de la categoría: Lo que he aprendido

Lo que he aprendido: de SVG a TikZ con Inkscape

Amo Inskape. A pesar de estar en la versión 0.92, funciona maravillosamente, tiene montones de funcionalidades y, por si fuera poco, al más puro estilo software libre, se pueden escribir extensiones que hagan lo que nos apetezca.

Como yo no soy tan hábil como para escribir mis propias extensiones, me aprovecho de la sabiduría de la comunidad para hacer cosas como la que os venía a contar hoy: exportar un svg a TikZ.

Un poco de contexto

Inkscape carrula con imágenes en formato SVG (Scalable Vector Graphics), un formato de gráficos vectoriales definidos en un XML. Que son texto plano, vamos, y podemos abrir una imagen SVG con un editor de texto y ver su definición. La ventaja de usar gráficos vectoriales es que no perdemos calidad al ampliar las imágenes porque los elementos se definen en función de su posición en el dibujo, si el dibujo es más grande las posiciones cambian y los elementos se dibujan de acuerdo a ellas.

TikZ, por su parte, es un paquete de TeX para crear gráficos vectoriales programáticamente1 (qué palabra más maravillosa). De esta manera, los gráficos nos quedan perfectamente integrados en nuestro documento y podemos cambiar el tipo de fuente, tamaño o lo que sea de ambas cosas simultáneamente.

Bien, la cuestión aquí es que usar TikZ es tremendamente difícil y usar Inkscape tremendamente sencillo y ¡LaTeX no tiene soporte para SVG! En general, suelo exportar las imágenes de Inkscape como EPS, manteniendo el texto como LaTeX, pero hay veces que me resulta más práctico tener la definición de la imagen en el propio archivo con el texto.

La extensión svg2tikz

Para esos casos en los que por lo que sea queremos que nuestra imagen esté definida en TikZ, tenemos la extensión de Inkscape svg2tikz. Son tres archivos, dos de extensión de Inkscape y uno de Python que están en la carpeta extensions del repositorio y que tenemos que copiar en la carpeta de extensiones de Inkscape, tal y como indican en las instrucciones. Será o bien /usr/share/inkscape/extensions o ~/.config/inkscape/extensions/2.

Cuidado aquí que yo guardé los archivos directamente del GitHub con Guardar destino como… y luego no me reconocía la extensión. Descargando todo el repo y copiando los archivos correspondientes en cambio sí que me funcionó.

Ahora en el menú Extensiones > Exportar tendremos una nueva opción Export to TikZ path que nos permitirá guardar la imagen como TikZ. Vigilad la ruta porque por defecto no la guarda en la carpeta donde estaba la imagen original sino en la carpeta de las extensiones.

Y ale, ya podemos chulear diciendo que nosotros hacemos todo directamente en LaTeX como hacían los profesores de mecánica de mi universidad que calculaban la base y la ruleta con la fórmula y luego decían que se veían directamente en el dibujo.

Referencias

Chapter 22. Extensions en el manual de Inkscape

svg2tikz

Galería de extensiones de Inkscape

inkscapeMadeEasy


Suena:


  1. Lo usé de tapadillo para pintar estructuras y que me quedaran cuquis lo apuntes, no sé si recordáis. 
  2. Como uno de mis propósitos del 2019 es ir abandonando las tecnologías privativas, voy a intentar hablar lo menos posible de Windows. 

Lo que he aprendido: hourglassing

No sé si os lo he contado alguna vez, pero mi investigación actual implica escribir código de elementos finitos. No, no son elementos muy delgaditos, sino finitos en contraposición a infinitos, es decir que su tamaño es limitado.

Otro día si eso me detengo un poco más en la explicación de los cacharros, que son movidas locas matemáticas, hoy lo que quería es hablar de un problema que se da a menudo en los cálculos de elementos finitos: el hourglassing. Es importante saber a qué se debe y cómo solucionarlo si queremos que nuestra simulación se parezca a la realidad.

Hourglassing: modos sin energía

El hourglassing (relojdearenización) consiste en que aparezcan modos de energía cero al utilizar integración reducida en un elemento. O, dicho en otras palabras, que el elemento sea más blando de la cuenta y se deforme sin que nadie lo toque. De esta manera, cuando calculemos los modos de nuestra estructura, nos aparecerán un montón de modos con un frecuencia cercana a cero. Dependiendo del orden del elemento, habrá más o menos modos de este tipo porque podrá deformarse de diferentes maneras.

Hay que tener cuidado aquí y no pensar que siempre que tenemos muchos modos de frecuencia cero hay hourglassing, ya que puede ser que hayamos metido la pata definiendo las condiciones de contorno, tengamos alguna pieza suelta y sean modos de sólido rígido perfectamente correctos.

Integrando en el elemento

Cuando digo integración reducida, me refiero a usar menos puntos de integración de Gauss de los que nos darían el resultado exacto. Por ejemplo, para un hexaedro lineal, necesitamos dos puntos de integración en cada dirección para obtener una matriz de rigidez exacta, a esto lo llamamos integración completa. Si decidimos usar solo un punto de integración en cada dirección, nuestro cálculo será menos exacto pero más rápido y corremos el riesgo de tener hourglassing. Ocurre lo mismo en el caso cuadrático, con tres puntos en cada dirección tendríamos integración completa y con dos integración reducida.

Podría pensarse que si la integración reducida da esos problemas es absurdo usarla, pero no es así, ya que nos evita un problema que padecen los elementos lineales: el shear locking. Esta otra palabra raruna hace referencia a la aparición de una fuerza cortante de mentira al intentar deformar a flexión un cuadrado o cubo lineal. Como es lineal, sus aristas no pueden curvarse, se comporta como un elemento mucho más rígido y nos da resultados absurdos. Aquí el bailoteo que permite la integración reducida juega a nuestro favor.

Visualizando el hourglassing

Aprovechando mi propio programita, os voy a enseñar la pinta que tienen los resultados de una simulación cuando ha habido hourglassing. En este ejemplo concreto es fácil de ver, ya que la superficie queda como arrugada. También se entiende mejor de dónde viene el nombre.

hourglass

Eso sí, cuidado, ya que no siempre es tan evidente. Lo infalible es calcular los modos: aparecerán los modos reales de la estructura y los espurios de energía cero que no tienen contrapartida física.

Conclusiones

Visto lo visto, os dejo con unos consejos para cuando simuléis, igualmente válidos para los que estéis escribiendo código desde cero y los que uséis software comercial:

  • Si salen modos con valores raros, después de verificar que las condiciones de contorno están bien, yo miraría qué tipo de integración estamos usando.

  • Si cambiando a integración completa los valores siguen sin tener sentido, yo comprobaría el orden del elemento. Si estamos usando elementos lineales con integración completa en un caso de flexión, son demasiado rígidos y van a dar valores de frecuencias naturales muy por encima de los correctos o deformaciones mucho más pequeñas (¡shear locking!). Una opción es pasarse a elementos cuadráticos, aumentando considerablemente el tiempo de cálculo, o considerar los elementos de modos incompatibles, que evitan el shear locking añadiendo restricciones adicionales.

  • Siempre verificaría que los problemas no se deben al tamaño de la malla, no cuesta (casi) nada reducir el mallado y nos quita muchos quebraderos de cabeza. Los elementos chiquiticos son más rígidos y menos propensos al hourglassing.

  • Si todo lo anterior falla, si estuviera en un programa comercial hurgaría en los controles del hourglass o metería amortiguamiento (si estoy en un caso que aplique). Solo como último recurso, todas las mierdas artificiales que le metamos a nuestra simulación tienen tendencia a explotarnos en la cara.

Espero que os resulte útil y poder algún día compartir el código para que juguéis.

Referencias

Why worry about hourglassing in explicit dynamics. Part 1

Why worry about hourglassing in explicit dynamics. Part 2

FEM: Rank deficiency and hourglassing

What is shear locking in FEA?

Lo que he aprendido: ejemplos con código en LaTeX

Estaba currelando en mis traspas para el curso de LaTeX para los doctorandos de mi uni y me pareció buena idea mostrar junto el LaTeX y el resultado para las tablas, ecuaciones y demás mandanga. Sabía que el paquete tcolorbox permite hacer esas historias, pero me parecía demasiado para la poca cosa que quería hacer. Así que hurgué un poco y me encontré con showexpl, un paquetillo que sirve precisamente para hacer lo que yo quería.

No os voy a contar gran cosa porque no deja de ser un añadido a listings que nos da el entorno LTXexample para escribir el código en un lado y el resultado en el otro. De hecho, en las opciones del entorno se puede meter la configuración de los listados y los pies funcionan exactamente igual que los de lstlisting.

Os dejo un ejemplo con los colores de la plantilla de beamer que estoy usando para que juguéis, es bastante claro lo que hace todo:

\documentclass[a4paper]{article}

\usepackage{xcolor}
\usepackage{amsmath}

\usepackage[final]{showexpl}

% Definición de colores
\definecolor{colororange}{HTML}{E65100} % orange
\definecolor{colordgray}{HTML}{795548}  % dark gray for note
\definecolor{colorhgray}{HTML}{212121}  % heavy dark gray for normal text
\definecolor{colorgreen}{HTML}{009688}  % green
\definecolor{colorlgray}{HTML}{FAFAFA}  % background light gray
\definecolor{colorblue}{HTML}{0277BB}   % blue

% Estilo de código
\lstset{
        tabsize=2, 
        backgroundcolor=\color{colorlgray},
        captionpos=b, 
        basicstyle=\small\ttfamily,
        columns=fixed, 
        extendedchars=true, 
        breaklines=true,
        prebreak = \raisebox{0ex}[0ex][0ex]{\ensuremath{\hookleftarrow}},
        showtabs=false, 
        showspaces=false,
        keywordstyle=\bfseries\color{colororange},
        commentstyle=\color{colorgreen}, 
        stringstyle=\color{colorblue},
        language={[LaTeX]TeX},
        texcsstyle=*\color{colororange}
}

\begin{document}

% Ejemplo de código
\begin{LTXexample}[pos=r,wide,width=.7,rframe={}]
% - pos: lugar donde saldrá el resultado
% - width: anchura relativa del resultado
% - rframe={}: para que no me ponga un cuadrado alrededor
\begin{equation}
  A = \pi\times R^2
\end{equation}
\end{LTXexample}

\end{document}

que nos crea esto:

A la izquierda el código de LaTeX necesario para escribir una ecuación. A la derecha la ecuación.

Ale, cambiadme parámetros y cosas y me decís qué os sale.

Lo que he aprendido: de Matlab a Julia

Hoy vengo a contaros una historia que se inició cuando volví de las vacaciones y descubrí que necesitaba reactivar la licencia de Matlab. Esto llevó a un lío en el que el servicio informático me desinstaló el programa y me instaló una versión hipervieja, en mi departamento me contaron que no, que no tenía que haber contactado con el servicio informático que la licencia la llevaba otra persona, tuve que desinstalar el Matlab viejo, reinstalar el nuevo y activar la licencia como decía esa persona. En definitiva, dos días de trabajo perdidos.

Si ya detestaba Matlab, esto me llevó a odiarlo aun más y buscar una alternativa libre y viable lo antes posible. Ya había andado con GNU Octave antes, pero no quería pasarme a algo que fuera siguiendo a Matlab (y a su sintaxis incoherente, sus funciones cambiantes de versión a versión y demás mierdas variadas que me sacn de quicio). Así que me puse a investigar lo que según el Fediverso es la mejor alternativa libre: el novedosísimo lenguaje de programación Julia cuya versión 1.0 es de agosto de este mismo año.

Sobre el lenguaje

Dicen sus creadores que cuando diseñaron Julia querían un lenguaje abierto (lleva licencia MIT porque se creó en el MIT), tan rápido como C, tan dinámico como Ruby, homoicónico como Lisp, con una notación matemática tan sencilla como la de Matlab, tan generalista como Python, tan útil para cálculo estadístico como R, con tantas capacidades para trabajar con texto como Perl y tan bueno juntando programas entre sí como un intérprete de comandos. Algo sencillo pero que satisficiera a los expertos, que fuera interpretado y compilado. Poquita cosa, vamos.

Lo presentaron en un paper, cosa que a mí ya me gana, y con unos números bastante espectaculares en el benchmark en el que pulen a Matlab en tiempos de computación y se acercan peligrosamente a C++, que toman como referencia. Del lenguaje poco más os puedo decir porque no entiendo nada más, solo que está escrito en una mezcla de C, C++ y Scheme, cosa que lo hace cuanto menos exótico.

En fin, centrándonos en cómo ponernos ya mismo a programar en Julia, no tenemos más que ir a su página de descargas y elegir nuestro método favorito. Yo, de momento, me he instalado el cacharro en el Windows del trabajo.

Diferencias con Matlab

En mi caso concreto que vengo de Matlab, lo que me interesa saber es cuánto de mi conocimiento actual puedo reciclar y la buena nueva es que ¡mucho! La sintaxis de ambos lenguajes es practicamente igual con algunas pequeñas diferencias a las que me he acostumbrado rápido, por ejemplo:

  • Los paréntesis solo se usan para pasar argumentos a las funciones, para elegir elementos de una matriz se usan los corchetes
  • No se copian las matrices, cuando hacemos A = B si cambiamos B cambia A
  • Para aplicar cierta operación elemento a elemento tenemos el operador splat ... o la construcción f.(x)
  • No se agrandan las matrices automáticamente, lo que nos obliga a prealquilar memoria
  • Se pueden tener todas las funciones que se quiera en un mismo archivo

Para los primeros días de locura tenemos hasta un traductor Matlab – Julia que nos ayuda a detectar estas pequeñas diferencias que son difíciles de ver.

Lo que más he notado es que Julia nos obliga a tener buenas prácticas a la hora de programar. Os pongo un ejemplo que me encanta: cuando calculas los autovalores con un método iterativo como Lanczos, te da cada vez cosas diferentes porque empieza con un vector aleatorio si tú, so vago matlabil, no le dices de dónde debe empezar, cosa que deberías hacer porque es tu matriz.

Otro tema chulo y diferente a Matlab es que Julia tiene un gestor de paquetes y debemos decirle explícitamente qué paquetes estamos usando en nuestro script o función con using PAQUETE. De momento he usado Polynomials, Plots, un paquete para cálculo simbólico y otro para matrices dispersas, todos ellos han funcionado muy bien y solo les puedo poner como pega que a veces la documentación es liosa, supongo que porque acaba de salir la versión 1.0. Se pueden además usar paquetes de Python, R o Matlab lo que hace que siempre tengamos recursos para nuestro caso de uso.

Mi entorno juliano actual

Ahora mismo estoy usando la consolilla de Julia (el Julia REPL) en Atom con la extensión Juno. Fui muy reacia a instalar Atom porque consume muchísimos recursos, pero ahora mismo en Windows es la mejor opción. Se instala en nada, tienes todo centralizado y funciona correctamente.

Y tiene un nivel de cuquismo elevado como podéis ver en la siguiente captura de pantalla:

Captura de pantalla de mi entorno para usar Julia en Atom

¿Por qué no desde Emacs?

Creedme cuando digo que lo intenté pero fue imposible. La explicación es sencilla: Windows. La mitad de los paquetes que permiten usar Emacs como IDE de Julia no carrulan en Windows. En teoría con instalar julia-mode y julia-repl o julia-shell-mode debería ser suficiente, pero después de horas infructuosas me rendí. Aun así, mereció la pena la investigación porque descubrí ESS, un modo de Emacs oara hacer análisis estadísticos al que espero volver algún día.

Conclusión

Como conclusión os digo que aunque lo he tocado poquito, Julia me encanta, hace lo que yo quiero, es natural programar en él y hay todo tipo de paquetes. Intentaré ir migrando todo mi código de Matlab a Julia e iré compartiendo con vosotros truquillos que vaya descubriendo.

¡En algún momento habrá código que Julia para calcular con elementos finitos gracias a mí! Hasta ese día podéis ir leyendo los enlaces que os dejo en las referencias 🙂

Referencias

Julia: A Fast Dynamic Language for Technical Computing

Julia for Matlab users

Crítica a Matlab en Neuroplausible

How I self-radicalised (an open science manifesto, part 1)

Documentación de Julia

A Deep Introduction to Julia for Data Science and Scientific Computing

Comparación entre Python, Matlab y Julia

Diferencias entre Matlab y Julia

Podcast con los creadores del lenguaje

Chuletilla

Think Julia

Julia Express (pdf)

Lo que he aprendido: dibujar cubos en Matlab

Como ando jugando con elementos finitos en el curro, me vi en la circunstancia de tener que pintar un modelo mallado con cubos en Matlab. Esto, que parece una tontería a simple vista es una locura total porque el buen Matlab hace solo superficies. Esto significa que para pintar un mísero cubo hay que pintar 6 superficies, cosa que tampoco es especialmente sencilla.

Después de probar varias alternativas, me he quedado con la función fill3 que pinta polígonos en tres dimensiones. Hay que meterle los datos de manera completamente antiintuitiva (es Matlab, recordemos), pero consigue el objetivo de crear cubetes.

Así que escribí esta función con la idea de usarla de base para pintar todo el modelo:

function cubo(coord)
% Pintar un cubos a partir de las coordenadas definidas según el siguiente
% orden para cada cubo:
%
%    7-------8
%   /|      /|
%  / |     / |
% 5--|----6  |
% |  3----|--4
% | /     | /
% 1-------2

for k = 1:length(coord)/8

  X = coord(8*(k-1)+1:8*(k-1)+8,1);
  Y = coord(8*(k-1)+1:8*(k-1)+8,2);
  Z = coord(8*(k-1)+1:8*(k-1)+8,3);

  % Los puntos de cada cara se ordenan según el sentido antihorario
  caras = [1 2 4 3; 5 6 8 7; 1 3 7 5; 2 4 8 6; 1 2 6 5; 3 4 8 7];

  % size(X) = [4 6]
  % - cada columna hace referencia a los puntos de un plano
  % - hay 4 elementos en cada columna que se refieren a las coordenadas x de
  % cada punto del plano

  X = [X(caras(1,:)) X(caras(2,:)) X(caras(3,:)) X(caras(4,:)) X(caras(5,:)) X(caras(6,:))];
  Y = [Y(caras(1,:)) Y(caras(2,:)) Y(caras(3,:)) Y(caras(4,:)) Y(caras(5,:)) Y(caras(6,:))];
  Z = [Z(caras(1,:)) Z(caras(2,:)) Z(caras(3,:)) Z(caras(4,:)) Z(caras(5,:)) Z(caras(6,:))];

  alpha = 0.6; % transparencia de la cara
  colour ='blue'; % color de la cara

  fill3(X,Y,Z,colour,'FaceAlpha',alpha); % dibujar los cubos
  axis equal

  hold on

end

end

Se usa de la siguiente manera:

P = [1 0 0; 1 1 0; 0 0 0; 0 1 0; 1 0 1; 1 1 1; 0 0 1; 0 1 1];

cubo(P)

Y nos pinta este simpático cubo:

cubo

Seguir leyendo

Lo que he aprendido: leer matrices MMF en Julia

Siguiendo con mi aprendizaje de Julia como sustituto a Matlab, he migrado una función que tenía para leer las matrices de masa y rigidez que exporto de Ansys para luego hacer movidas locas. Otro día os cuento el por qué de la transición (ejem, licencia, ejem) y cómo podéis uniros a mí en el lado oscuro.

Ale, al tema. Primero os voy a contar cómo exportar de Ansys en un par de párrafos, os doy un poco de contexto sobre escribir matrices en archivos de texto y luego vamos al grano.

Exportar matrices de Ansys

Para obtener las matrices de masa y rigidez de Ansys en primer lugar
exportamos el archivo input del modelo con: Tools > Write input file. Luego, a ese input le añadimos antes de la línea /wb,file,end este trocico:

! Gets Stiffness Matrix
! D for double precision real number, C for complex
*SMAT, MatK, D, import, full, file.full, stiff
*export, MatK, mmf, matkMMF.txt ! Exports Stiffness as MMF format

! Gets Mass Matrix
*SMAT, MatM, D, import, full, file.full, mass
*export, MatM, mmf, matmMMF.txt ! Exports Stiffness as MMF format

Seguir leyendo →

Lo que he aprendido: integración numérica con Gauss

Esto no es técnicamente algo que haya aprendido ahora, llevo años integrando numéricamente como manda Gauss, pero creo que siempre lo había hecho o bien a mano o bien buscando los puntos de integración y los pesos en las tablas. La novedad aquí es que uso los polinomios de Legendre para calcular cuáles son los puntos de integración y los pesos correspondientes para cualquier número de puntos.

Recordemos que en el método de Gauss:

  • Los puntos de integración corresponden a las raíces del polinomio de Legendre de orden n
  • Los pesos se obtienen a partir del valor de la derivada en los puntos de integración
  • El método se define para el intervalo [-1,1], para cualquier otro intervalo tenemos que hacer un cambio de variable y ajustar el valor final de la integral

Con mis reducidas habilidades julianas, me queda una función con esta pinta (está in inglis porque la migré de Matlab de otra función que es parte del código de elementos finitos que estoy fabricando en el curro):

# Pkg.add("Polynomials")
using Polynomials

"""
legendre(x)

Create Legendre polynomial of degree `x`
"""
legendre(x) = 1//(2^x*factorial(x))*polyder(Poly([-1, 0, 1])^x,x)

"""
gaussGeneral(fun, points, interval)

Compute numerical integral of function `fun` using number
of `points` in `interval`
"""
function gaussGeneral(fun, points, interval)

a = interval[1]
b = interval[2]

# Compute roots of Legendre polynomial of order 'points'
p = roots(legendre(points))
d = polyder(legendre(points))

# Weights
w = 2*ones(length(p),1)./((ones(length(p),1)-p.^2).*d(p).^2)

# Variable change
p = (a+b)/2*ones(length(p),1) + (b-a)/2*p

# Sum of weighted values
return (b-a)/2*w'*fun(p)

end

Lo chachi de definir así la integral numérica es que me permite calcular la integral de la función con diferente número de puntos y comparar con el valor analítico. Yo soy ingeniera y no me creo las típicas frases la cuadratura de Gauss es un método para integrar de manera exacta polinimios hasta de orden 2n-1 donde n es el número de puntos utilizados si no lo pruebo.

Y eso era. Dejo una pregunta para los sabios: ¿me podéis recomendar una referencia en la que se explique de manera entendible por qué las raíces del polinomio de Legendre son los puntos óptimos de integración? Creo que nunca lo llegué a comprender.


Suena ahora: