The General Curl
The Theory
The image you posted above is trace visualization for what is known as a divergence free vector field. To understand what that means, consider particles that move along the lines you see above (the field represents their instantaneous velocity), if the field is indeed divergence free, those particles will never collide. That's what gives those visualizations their beauty, the lines never intersect.
Mathematically, a vector field $F$ is divergence free if $\nabla \cdot F = 0$. It is hard to compute a random vector field that is divergence free proceduraly, so we use some math to aid us. In vector calculus, the Curl of any vector field is divergence free. The curl is an operator that operates on a vector field and returns another vector field that represents the infinitesimal (microscopic) rotation of the input vector field. The infinitesimal rotation of a vector field can be understood using this intuitive interpretation from Wikipedia1:
Suppose the vector field describes the velocity field of a fluid flow (such as a large tank of liquid or gas) and a small ball is located within the fluid or gas (the center of the ball being fixed at a certain point). If the ball has a rough surface, the fluid flowing past it will make it rotate. The rotation axis (oriented according to the right hand rule) points in the direction of the curl of the field at the center of the ball, and the angular speed of the rotation is half the magnitude of the curl at this point.
We are not interested in what the curl represents, we are only interested in the fact that the curl of any vector field is divergence free. For our artistic use, the input of the curl operator $F$ can be a simplex or perlin noise-based vector field, that is, $F = (F_x, F_y, F_z)$ where $F_x, F_y, F_z$ are different perlin or simplex noise functions.
Let us now look into the curl operator that will take $F$ as an input and returns a divergence free vector field, the curl of $F$ denoted by $\nabla \times F$ is equal to:
$$
\nabla \times F =
\left( \frac{\partial F_z}{\partial y} - \frac{\partial F_y}{\partial z} \right)\mathbf{i} +
\left( \frac{\partial F_x}{\partial z} - \frac{\partial F_z}{\partial x} \right)\mathbf{j} +
\left( \frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} \right)\mathbf{k}
$$
Where $\mathbf{i}$,$\mathbf{j}$ and $\mathbf{k}$ are the fundamental vectors of the rectangular space. Don't worry if you don't understand those equations, we shall understand them soon. The notation $\frac{\partial F_z}{\partial y}$ represents the partial derivative of $F_z$ with respect to the $y$ axis, similarly, $\frac{\partial F_y}{\partial z}$ represents the partial derivative of $F_y$ with respect to the $z$ axis, and so on. A partial derivative measures how much a function changes along an axis. Take a brief look at derivatives before we continue, you don't have to study it, just recognize what it represents. To compute the partial derivative, we use what is known as the central finite difference method or Symmetric Derivative which states that the partial derivative is equal to:
$$
\frac{f(x+h) - f(x-h)}{2h}
$$
Where $h$ is some arbitrary small number and $f$ is the function in question. The same equation applies for the $y$ and $z$ axis. For instance to compute $\frac{\partial F_z}{\partial y}$, we add $h$ (a small number) to the $y$ component of the the vector, evaluate the $F_z$ noise at that vector, subtract $h$ from the $y$ component of the vector, evaluate the noise and take the difference between the two evaluated noise values, and finally divide them by $2h$.
To understand how the central finite difference method work, consider this intuitive interpretation:
Suppose one is standing on a mountain blindfolded, one was asked to report the steepness (slope, partial derivative) of the point one is standing at, one might move a step to the right and identify whether one ascended or descended and how much did one ascend or descend by, another might move a step to the left and identify whether one ascended or descended and how much did one ascend or descend by. Now consider the situation where one is at the very top of the mountain, at that point, the slope is zero, one is not ascending or descending, but both persons will report that they are descending, because if one move to the left or the right one will be descending. We conclude that to get an accurate result, we have to move a step to the right, observe, move a step to the left, observe and then determine if you are ascending or descending. In our case, if one move both to the left and to the right, one will understand that one is at an even surface (at the top). So when we say $f(x+h)$ we simply mean that we observe after we move a step $h$ to the right and $f(x-h)$ means to the left.
Once we compute the curl of the vector field, we can do what is known as advection using Euler's integration. When treating the curl vector field as the instantaneous velocity of some particles in space, Euler's integration approximates the location of the particles after one second by adding the curl vector field to the original particles location. So to get the whole trace of a particle, we evaluate the curl, add it to the location of the particle, evaluate the curl at the new location, add it to the new location, and so on.
Implementation
Simple Curl Implementation
The first step is to make a group that generates the vector field by initializing its components with different simplex noises $(F_x, F_y, F_z)$:
The next step is to compute the partial derivatives. This can be done using the formula for the symmetric derivative:
Take your time to understand this node tree. It can be hard to grasp at first. Having all the partial derivatives, we can now compute the curl easily using the equation above:
And by advecting some initial vectors along the curl vector field using Euler's integration and creating splines from the output points:
We get something like this:
Optimized Curl Implementation
While the previous implementation works, it is very slow because it requires a lot of noise evaluations and noise is expensive to compute. Thankfully, the noise functions in Animation Nodes v2.1 is implemented using SIMD instructions, which means it can somewhat be executed in parallel for CPUs that support Advanced Vector Extensions (AVX) which is pretty much every modern CPU. This result in a speed up of up to 600x the original implementation. SIMD requires the input data to occupy a contiguous memory block. So to utilize SIMD instructions and optimize this setup, we should combine all our data into a single big vector list, evaluate the noise at it, segment the output, and process it. The following implementation is hard to understand, so it is better to do it on your own, my implementation looks like this:
Euler's integration loop stays the same:
However, notice that we are appending a list of vectors and not a vector, so the output will be of length $n \cdot m$, where $n$ is the number of iterations and $m$ is the number of input initial vectors. We have to do another segmentation before creating splines from those points, so the spline loop will include segmentation and spline creation:
And what we get is a fully functional, fast, curl trace generator:
Blend file for study and practice:
Surface Curl
Now that we have an understanding of what the curl is, we can go ahead and finally answer your question.
Theory
Let us define the problem as follows:
We want the output vector field of the curl operator to always be tangent to the surface of our meshes, that is, the curl should always be perpendicular to the vector field that represents our mesh surface normals.
If our initial points were on the surface of the mesh and the above condition is satisfied, then we will get the result we are looking for. The above condition is satisfied only if $(\nabla \times F) \cdot \vec{N} = 0$, where $(\nabla \times F)$ is the curl and $\vec{N}$ is our normals vector field. The problem is now to compute a vector field $F$ so that the dot product between its curl and the normals vector field is equal to zero.
By some analysis we get, $F$ is equal to $\vec{N} \cdot p(x,y,z)$ where $p$ is a perlin or simplex noise function just like we used in the general case, the curl in that case becomes:
$$
\begin{aligned}
F &= p(x,y,z) \begin{bmatrix} N_x \\ N_y \\ N_z \end{bmatrix}\\
(\nabla \times F) &=
\left( \frac{\partial F_z}{\partial y} - \frac{\partial F_y}{\partial z},
\frac{\partial F_x}{\partial z} - \frac{\partial F_z}{\partial x},
\frac{\partial F_y}{\partial x} - \frac{\partial F_x}{\partial y} \right)\\
&= \left( N_z \frac{\partial p}{\partial y} - N_y \frac{\partial p}{\partial z},
N_x \frac{\partial p}{\partial z} - N_z \frac{\partial p}{\partial x},
N_y \frac{\partial p}{\partial x} - N_x \frac{\partial p}{\partial y} \right)
\end{aligned}
$$
And the proof:
$$
\begin{aligned}
(\nabla \times F) \cdot \vec{N} =
&N_z N_x \frac{\partial p}{\partial y} - N_y N_x \frac{\partial p}{\partial z} +\\
&N_x N_y \frac{\partial p}{\partial z} - N_z N_y \frac{\partial p}{\partial x} +\\
&N_y N_z \frac{\partial p}{\partial x} - N_x N_z \frac{\partial p}{\partial y} = 0
\end{aligned}
$$
Notice that if the normal is $\vec{N}=(0,0,1)$, you only end up with the equation that Robert Bridson proposed in his paper "Curl-Noise for Procedural Fluid Flow", which is:
$$
\frac{\partial p}{\partial y}\mathbf{i} -\frac{\partial p}{\partial x}\mathbf{j}
$$
So now you know where my other answer came from. Since the curl have zero z component, we shall call this special case the 2D Case.
Implementation
The implementation is much easier because we only have to differentiate a single function. You should by now know how to implement that. Here is the non-trivial implementation of the surface curl; notice how there is an extra input called normals now:
Lets consider the simple example in which we use a sphere as our mesh, the normal of the sphere at any point on its surface is equal to that point normalized, so the normal is simply the point locations normalized. After the curl is computed and the the points advected, we can project the output onto the surface to make sure the splines always lie on the surface of the sphere since we are merely performing a finite approximation. We won't need this step for sufficiently small step sizes for Euler's integration which is in our case decided by the magnitude of the normal vector field or the amplitude of the noise function. To project any point onto the surface of a sphere of radius $r$ we normalize the vector and multiply by the radius, so the advection loop looks like this now:
The result are those magnificent and beautiful splines:
Blend file for study and practice:
Since the 2D case have two zero normal components, let's implement the 2D case using an optimized setup. We simply remove anything that got multiplied by zero to get:
This gives:
Blend file:
Ok, I have given you an example where the mesh is defined implicitly, but how do we use an actual mesh? Well, we can use a BVH tree:
We approximate the normal field by the Nearest Surface Point node Normal output, and we change the generator to be the location of the nearest surface point. Why? Because this is the projection of the point on the surface, which is a step we have to do to make sure we get accurate results as we discussed before. Make sure to reduce the step size, because BVH trees are not as accurate as the implicit definition of our sphere.
It should be noted that constructing BVH trees from objects directly won't apply modifiers, so such node tree can be used:
And finally, here is the result of both setups:
Inviscid Boundary Condition
An inviscid boundary condition required the dot product between the curl and the surface normals around the boundaries of objects to be zero, in other words, it requires the curl to be orthogonal to the normals of the object around its boundaries. Ensuring this condition in the surface curl case is easy, while it is somewhat challenging for the general 3D case, so I shall only introduce the surface curl.
Looking back at the equation for the 2D case:
$$
\frac{\partial p}{\partial y}\mathbf{i} -\frac{\partial p}{\partial x}\mathbf{j}
$$
We notice that the curl is actually orthogonal to the Gradient, the gradient vector being, in simple terms, a vector pointing in the direction of the greatest change. The gradient is computed by:
$$
\frac{\partial p}{\partial x}\mathbf{i} +\frac{\partial p}{\partial y}\mathbf{j}
$$
Can you see the similarities? So we can implement the boundary condition by making sure the noise field have a gradient tangent to the normals of the object around its boundaries, in other words, the noise field should increase the most when moving in the direction of the surface normals. We actually know a field that satisfied this condition everywhere, that is, its gradient is tangent to the surface normals. This field is the Signed Distance Field (SDF). An SDF is, for every point, the distance to the closest surface point, having a negative sign if it is inside the surface and positive if outside (However, the sign is redundant in this particular application). So the condition can be implemented by modulating the noise field based on the SDF of the surface using a ramp. As Robert implemented it, the noise field should be multiplied by the ramp:
$$
\begin{cases}
1 & r > 1 \\
\frac{15}{8}r - \frac{10}{8}r^3 + \frac{3}{8}r^5 & -1 \leq r \leq 1 \\
-1 & r < -1
\end{cases}
$$
Where $r$ is the SDF divided by a scalar defining the modulation width. Which we can implement as follows:
What we did is replace the noise node with the above subprogram and add a BVH tree as an input.
Using Older Animation Nodes Versions
The only node used in this answer and does not exist prior to Animation Nodes v2.1 is the Vector Noise node which provides us with the perlin or simplex noise functions.
To make the node tree work for older versions, you simply have to replace the Vector Noise node with the built-in blender noise functions that exists in the mathutils python module. To do so, we can use the expression node:
With this code:
[mathutils.noise.noise(x) for x in vectors]
Aside from noise.noise, there are other functions that give you more control, which you can find in the API.
File for study:
