maanantai 23. helmikuuta 2015

Matlab EPS-export fails on dotted lines

I continue my set of tips on drawing bearable figures using Matlab. Have you ever tried to increase the line width of a double dotted (with linespec ':') or dot-dash (with linespec '-.') line? Does it look like crap after exporting to *.eps?

There is a simple solution to the problem: open up the *.eps-file and search for

  1. "/DO" for ':'
  2. "/DD" for '-.'
The corresponding lines should be something like


  1. /DO { [0.5 dpi2point mul 4 dpi2point mul] 0 setdash } bdef
  2. /DD { [0.5 dpi2point mul 4 dpi2point mul 6 dpi2point mul 4


By increasing the first decimal number 0.5 up to e.g. 1 you should get much better results.

keskiviikko 10. joulukuuta 2014

Trying out the Schrödinger

Time-independent Schrödinger equation constitutes the following linear eigenvalue problem: find $u : \mathbb{R}^3 \rightarrow \mathbb{C}$ and $\lambda \in \mathbb{C}$ such that
$$-\frac{h^2}{4\pi \mu} \Delta u + V(\boldsymbol{x}) u = \lambda u,$$ where $h$ and $\mu$ are some constants and $V(\boldsymbol{x})$ is a given potential. For hydrogen atom we have
$$V(r) = -\frac{e^2}{4\pi \epsilon_0 r},$$ where $e$ and $\epsilon_0$ are some additional constants and $r$ is the distance from the nucleus. Writing all the constants as one constant
$$C = - \frac{\mu e^2}{\epsilon_0 h^2}$$ and rescaling the eigenvalue leads to
$$-\Delta u - C r^{-1} u = \widetilde{\lambda} u.$$ Let $\Omega$ be the unit ball $B(0,1)$ and consider the homogeneous Dirichlet boundary condition $u|_{\partial \Omega} = 0$. Weak formulation for this problem is: find $u \in H^1_0(\Omega)$ such that
$$(\nabla u, \nabla v)-C(r^{-1}u,v)=\widetilde{\lambda} (u,v)$$ holds for every $v \in H^1_0(\Omega)$.

Setting $C=1$ and solving this with piecewise linear elements on my Matlab FEM code gives eigenvectors whose midsections look as follows:

First 12 lowest eigenvectors of the Schrödinger equation with constants normalized to 1.

You may compare them to the solutions given at Wikipedia. Can you find them all?

lauantai 6. joulukuuta 2014

On ray marching

One way to define a surface in three-dimensions is through a function $f:\mathbb{R}^3\rightarrow\mathbb{R}$ which satisfies
$$f(\boldsymbol{x})=0,$$ for every point $\boldsymbol{x}\in \mathbb{R}^3$ that belongs to the surface. Thus, the surface $S$ can be defined as
$$S=\{\boldsymbol{x} \in \mathbb{R}^3\,|\,f(\boldsymbol{x})=0\}.$$ For one to visualize the surface $S$, the intersection point $\boldsymbol{p} \in \mathbb{R}^3$ of a ray---starting from the camera position $\boldsymbol{c} \in \mathbb{R}^3$ to some direction $\boldsymbol{d} \in \mathbb{R}^3$, $\|\boldsymbol{d}\|=1$---and the surface must be sought. After these intersection points are known for sufficiently many rays, a visualization can be constructed (some examples follow).

One way to find the intersection point $\boldsymbol{p}$ is through analytical calculations: the surface defined by $S$ is sometimes so simple that this calculation can be done by hand. An example that considers the intersection of a ray and a sphere can be found from Wikipedia. This is exactly what the traditional ray tracing is all about.

In a more general setting, no analytical relationship exists for finding the intersection point $\boldsymbol{p}$. In such a case, one must rely to numerical methods.

A proficient applied mathematician would immediately attempt to apply, e.g., Secant method to the problem: find $t \in \mathbb{R}$ such that
$$f(\boldsymbol{c}+t\boldsymbol{d}) \approx 0.$$ There is, however, a caveat. For a general $f$ the function
$$g(t):=f(\boldsymbol{c}+t\boldsymbol{d})$$ might have several roots and in Secant method (or similar) nothing guarantees that a root nearest to the camera is recovered. You might end up drawing the wrong side of the object of interest! Thus, the applied numerical method should be guaranteed to find a root nearest to the camera.

A somewhat trivial method to find the nearest (or approximately nearest) root is to fix a sufficiently small constant step size $\Delta t$ and then find the smallest $n \in \mathbb{N} \cup \{0\}$ such that
$$f(\boldsymbol{c}+n \Delta t \boldsymbol{d})\approx 0.$$ This is known as ray marching and is indeed a valid method for visualizing the surface $S$. The problem is that it's tremendously slow since one must try every possible $n$ starting from zero!

Fortunately, there exists a clever trick to make the computation faster by ensuring that $f$ defines a distance field of $S$. This means that at every point $\boldsymbol{x}$ in the space $\mathbb{R}^3$ the value $f(\boldsymbol{x})$ represents the minimum distance from $\boldsymbol{x}$ to the surface $S$. This sounds like a hard restriction to satisfy but for many fractals the distance fields (or approximate distance fields) have already been derived by enthusiasts.

The algorithm goes as follows:

 1. $\boldsymbol{p}:=\boldsymbol{c}$
 2. $\boldsymbol{p}:=\boldsymbol{p}+f(\boldsymbol{p})\boldsymbol{d}$
 3. If $f(\boldsymbol{p})$ is larger than some tolerance go back to step 2. Otherwise, break.

Output of the prescribed algorithm is the intersection point $\boldsymbol{p}$.

If you happen to own a sufficiently fast graphics processing unit and know the function $f$ for some fractal then this method can be implemented in real time using, e.g., OpenGL Shading Language. See the following Shadertoy entry I constructed as an example:

Visualizing a mandelbulb fractal. See https://www.shadertoy.com/view/MtX3zr.

Getting the coloring and other visuals right requires of course lots of experience (or trial-and-error) but this is the basic idea for generating such animations. I tried to add some comments to the shader code and attempted to use same symbols as in this post.

tiistai 4. marraskuuta 2014

P.G. Ciarlet on the words "linear triangle" and "biquadratic rectangle"

Since the readers of this article are perfectly aware that a triangle may be isoceles, but certainly not genuinely "linear", or that a rectangle may be square, but certainly not genuinely "biquadratic", I have blithely committed these serious abus de langage, which clearly convey more information than the names that I had originally chosen.

P.G. Ciarlet, Handbook of Numerical Analysis, Part 1

maanantai 18. elokuuta 2014

Exporting to vectorized PDF images in MATLAB

Keywords: painters, pdf, trisurf, white lines, patch

Ever tried to export nice looking, smoothly shaded pdf figures from MATLAB and ended up with something like this?
Would you prefer something like this?
If the answer to both questions is yes, then read on.

A plenty of forum posts in Mathworks' website propose different kinds of magic tricks to solve this problem. However, none of the tricks seemed to work for my case where I was drawing a surface with trisurf and exporting it to vectorized PDF with the following code.
trisurf(mesh.t',mesh.p(1,:)',mesh.p(2,:)',...
                  ufun(mesh.p(1,:),mesh.p(1,:))');
shading interp;
print -painters -dpdf -r600 vel.pdf
This always ended up with the annoying white lines as demonstrated in the first figure. Finally I realized that trisurf uses underneath the patch command to actually draw the triangles. I came up with a simple fix of setting up interpolated edge colors to every triangle drawn by the patch command as follows:
H = trisurf(mesh.t',mesh.p(1,:)',mesh.p(2,:)',...
                  ufun(mesh.p(1,:),mesh.p(1,:))');
set(H,'EdgeColor','interp');
shading interp;
print -painters -dpdf -r600 vel.pdf
This gives the beatiful plot in the second figure. For cropping the whitespace around the pdf I propose the usage of the linux command line tool pdfcrop.

lauantai 17. toukokuuta 2014

3D progress

Today I added preliminary support for 3D bilinear and linear forms using tetrahedral P1 elements. Unfortunately, the visualization is based on MATLAB using MATLAB.jl Julia package (https://github.com/lindahua/MATLAB.jl) and the mesh generation is also based on a MATLAB script package iso2mesh (http://iso2mesh.sourceforge.net/cgi-bin/index.cgi) which is again based on tetgen (http://wias-berlin.de/software/tetgen/).

Thus, I haven't yet figured out how I should package the code properly for distribution. Meanwhile, here is a screenshot of solving the Poisson equation
$$\Delta u = 1$$ on a unit cube $[0,1]^3$ with zero Dirichlet boundary conditions and the corresponding mesh:

Cross section at $z=0.5$.
The mesh. Color indicates simply the $z$-coordinate.

sunnuntai 13. huhtikuuta 2014

Improved assembly

Yesterday I learned how to do variable argument macro's in Julia and based on my new knowledge implemented the bilinear assembly (and linear assembly) with material parameters using them. Thus, now there is only a single macro for bilinear assembly and bilinear assembly with a material parameter (or some other node-wise defined parameter). Next goals for the Julia-code:

  • Assembly of the boundary terms
  • Faster mesh refinement
  • Simple geometry representations for 2D-domains
  • 3D-linear tetrahedral element
  • ...
  • and in a distant future
  • ...
  • Geometry tools
  • Additional elements: p-refinement, vectorial elements (i.e. Raviart-Thomas)
  • Isoparametric mappings
  • Simple geometry representations for 3D-domains
Some of these might be too much to realize in practice but there is always hope.