today we're gonna start working on a big one
Questan's equation
and in order to attack this and give you a good
intuition I want to approach this in steps
and the way we're gonna talk about this
and think about intuitions
and what one do is based on Matt Ferrero's brilliant
block post called Persons Equation
is the most powerful tool
not yet in your toolbox
link in the description
I highly recommend reading this article
and in general
Matt Ferraro's Block is a treasure trove of physics
math and engineering
so before we get to Persons Equation
let's start with something a bit simpler
which is called LA Plaza's Equation
named after plc Mon LA plus
so LA Plaza's equation usually is written as this here
what does that mean
the first thing you'll notice is this
inverted triangle
that is squared at the beginning
and this is quite confusing
admittedly if you've never seen it
so this triangle here is called nabla or del
and this formula says that
nabla squared or del squared
times a function should be equal to 0
let's go with this
so this equation sometimes is also written like this
here
this is upright triangle being the LaPlace operator
so let's talk about
this Dell square and the LaPlace operator first
what do they mean or how do they come about
so this laplacean
the LaPlace operator
can be written as nabla times nabla or a del times del
or in brief del squared so nothing new here
just to let you know
why this
the flashing operator can also be written as this
del squared but what does it actually mean
this del times del or nabla times nabla
it actually is a functions second derivative
so let's go over this
the first derivative of a function so a single
now BLA is a functions gradient
and taking the derivative of a functions gradient by
again multiplying it with
bell will give us the divergence
of that gradient
so what we're saying in this formula
is that we want the divergence
of the gradient of our formula to be zero everywhere
again this
expression here means the diversion of the gradient
so what does that mean
finding a function whose divergence
of its gradient is zero everywhere
so let's talk about what a gradient is and what
divergence is
let's talk about the gradient first
what I have here is this curved wavy surface
and the surface
could be viewed as a function in three dimensions
so I feed in an X and a Z coordinate
and it spits out a y coordinate
which we plotted here on this grid
now the gradient of a functions point is the
direction of the steepest ascent
for that point
to give you an analogy
imagine hiking in a very foggy terrain
and you want to find the peak of a mountain
but you cannot see far just say half a meter or so
so the best strategy you could come up with
is just look around in your direct vicinity
and go into the direction
where the terrain rises most deeply
and that is the gradient
the direction of the steepest ascent
we can measure this using a measure node in Houdini set
to measure the gradient
in this case of our
why component
resulting in this feel of arrows a vector field
and if we zoom in
you can see
that all those arrows point in the direction
of the steepest ascent
for each of those points here making up our terrain
but in this case they all point to this peak
to this ridge here
when we look at the valleys in here
you can see that all those arrows point away
from those points in the valley
just trying to move uphill
and also you can see that the measure node here
colored are arrows
indicating the gradient
and you can see that those arrows get blue
in the areas where there are peaks or valleys
and in mathematical notation or mathematical language
those peaks or valleys are called
local maxima if we have peaks
or local minima if we have valleys
so on peaks and valleys we got those blue arrows
and on smoothly
almost uniformly rising areas of our mountain
we get those red arrows
so in here where we have those red arrows
all the surrounding arrows for one point point
nearly in the same direction
while in the blue areas up here
for one single point
all the surrounding arrows point
into wildly different directions
same goes for those valleys in here
if we zoom in we can see that for example
for this point here
the surrounding arrows point in
wildly different directions
turns out Houdini wants to visualize and help us
by showing us
the divergence of that gradient
what does that mean so what does divergence
let's have a look at a point in a vector field
the vector field could be gradient
or it could be velocities of a gas
and in this example here
the velocities surrounding our point all point
away from it
this could be the case for a gas being heated
at this very point here
so the gas is expanding
influence simulations
this behavior sometimes called a source
so calculating the divergence
of the vector field at this point here results in a
positive number
meaning that the arrows
or the vectors of that vector field
point away from this point here
in contrast we could look at this configuration here
where all velocity surrounding that point
point towards it
in this case
we could think of it as gas being cool at this point
and thus contracting and influence
sometimes it's called a sink
and if we calculate
the divergent of this vector field here
at this point
we get a negative number
meaning that the vector field
points at this point here
contracting towards it
and finally in this case
the velocities surrounding our point
all point in the same direction
could think of gas flowing uniformly at this point
and if we calculate the divergence
of the vector feel at this point here
we get zero and that's exactly what's visualized here
in Houdini's measure node with those colors
so back to our equation what does all of this mean
well this expression here reads as
find me a function
where the divergence of the gradient is
zero at every point
again if the diversions of the gradient is non zero
such as in these blue areas here
it does mean our function has local minimum maximum
and that's what we definitely don't want
rather what we want is a function
that is as smooth as possible everywhere
except at its boundaries
that means its edges
to give you an idea of what I mean
or how such a function might look like
I prepared this grid here
quite high resolution on it
I prepared a boundary group so just these edges here
and I specified one single edge group here
this one onto which I write a high value
in this case one
and the rest of the grid has a value of 0
so if I plot those values along the wide erection
you can see that only this edge here has high values
while the rest of the grid
has a value of 0 on every point
and now
I want to fill the rest of the grid in such a way
that the divergence
of the gradient of the values that are right in here
is as close to zero as possible
that means that those values which I will
calculate in here will change as smoothly as possible
in this case
I can use Houdini's new attrib fill note to do that
and the resulting function looks like this
so it's this curved minimal surface that spans between
those boundary edges here
let's see how well we minimize
the diversion of the gradient
by just measuring the gradient again
and visualizing it in here
let me just switch off
the underlying color visualization here
and you can see that except
for some areas in the boundaries
we did a decent job and kind of succeeded at making
the direction of this functions gradient
point roughly in the same
direction for every point
and its surrounding points
how does the attribute fill do that
or how do we do that
how do we solve this condition
that we want to fill this area as smoothly as possible
what we do is we simply take one point for example
let's take this one here
where our function suddenly gets steep
and we calculate its value as being
the average of the surrounding points values
so we try and average out the distribution
of that value
for every grid point
and this in my opinion is where it gets interesting
because there are several algorithms out there
to do that making your life easier or harder and being
faster or slower
so let's talk about them one by one
so what I'll be doing here is implement what's called a
Jacobi solver or Jacobi iterations using VEX first
we'll get to open CIA later
but for readability
and not to have to deal with open CIA specifics
I'm gonna go VEX first
and compare the solver to another technique called
gossile iterations
or gossile so what I want to do is drop down a grid
which will specify my function space
so where I want to set my values make it 1 by 1 units
and just for now for testing make it rather low
rest by 32 by 32 columns
grid of that grid here I think
also maybe safe this
and then as in the demo I'm going to create two groups
call one boundary disabled based group
set this one to be on points and include by edges
and down here I want to check unshared edges
so I get all those boundaries in that group
also I'm just going to create one group
where I set the values to 1
so if I struck down another group
said this to be points as well
disable the base group and check
keep in bounding regions
and in here I want to do the same thing as
it previously
just select this row points here
and then using an attric create
just for the group I just created
I'm going to specify the value
which in usual mass notation is just
called f and my default
I want to set this value to 0 everywhere but to 1
on the group I just created
let's just visualize that
by going to the info button here and clicking on the f
we just created
and we can see in this blue
and red color the values of our function here
let's right click on the visualizers here
and edit that gradient and
set it to be infrared which I find a bit more readable
next for calculating the average of a point
being the average value of its surrounding points
I need some attribute storing
which are at points given neighbors
which I can easily do using a point wrangle
and in here creating an array call it NPTs
and in here I want to store the neighbors
of our currently active point
done like that
also I want to know how many points I found
so how many neighbors a given point has
and let's call this attribute
num and should be equal to the length of that array
that list and pts we just created
let's store both as a point attribute
first gonna be an integer list
let's call this one npts
and should be equal to the NPTs we just found
and also let's store that integer attribute
which we gonna call num
and should be equal to num
let's just have a look in geospread sheet
got the neighbors
that is two neighbors that is three neighbors
okay that's looking good back in here
let's maybe zoom in a bit on that one
so let's finally start calculating
each point's new value
just taking the average of the surrounding points
and setting that to be each point's new value
in again using a point wrangle
to add here and then here
the first thing I wanna specify is the average
so that's the new value of our point in the beginning
you want to electrolyze this to be 0
and then I wanna run over all the neighboring points
so I'm gonna use a for each loop
using N P t s
my variable into which we store the current point
out of our list of neighboring point
which I can look up using the atrib short hand
like this remember that's the N P t s
we stored in the previous wrangle here
and then for each of those points
I wanna look up there individual value
so let's call this one N F
as equal to the
points F value stored in the neighboring point
and all we have to do is add this to our average
like so and at the end
after we summed up all those values
from the adjacent points
making sure
to divide that value by the number of adjacent points
so I'm going to say average divided equals and again
that's short hand
for average equals average divided by
I want to cast this to float
and I want to look at the number
of adjacent points which I
previously stored in that attribute here
in the previous point triangle
like so and now all I have to do is set that points
new F value to be this average
well something happened
that value creeped forward one row of points
when we compare this to the previous state
something happened
and again let's visualize this
and maybe before we do that
just call those point wrangles
something more descriptive than point wrangle 1 and 2
so let's call this one in it neighbors
and let's just call it in it neighbors
this one here let's call it Kalk Laplashim
and let's use another pointer angle
which we're going to call
visualize a value of f
and in here
all I want to do is offset my current points
y position via the points current value
so we can now see that after
running this calcula plash in here
let's just disable lighting here
before it was this state now it's this state
and we can see two things
on the one hand
this definitely isn't the function we'd expect
on the other hand
this boundary here is being manipulated
taken into account for calculation
so let's first exclude the boundaries by running this
over everything but the boundaries
exclamation Mark boundary selects this group here
so now that we have this calculate laplash
and running only over the interior points
we can see that
well
some of those points close to that initial value here
get pulled up a bit
the secret to calculating this
as you might have guessed is just running it
over and over and over again
in our case I'm just gonna use a for loop with feedback
I'm gonna cut this
and just move the very first note up in here
and the second note of that for loop down here
like so
so by default this is set to run over 10 iterations
and you can see it's already getting smoother
but we need quite more iteration so this use 64 28
to 56 5 12 and 10 24
and after 10 24
if I double this again to 2
48 you can see there's not much different
so this thing needs about 1,000 iterations
to spit out a decent solution for this
very simple formula
and clearly this takes time
although we might be
able to speed up our for loop using compiled locks
but still a lot of iterations
the upside of this whole algorithm is that in here
where we calculate each points new value
by just taking the average of its surrounding points
can be run in parallel
so the execution order meaning
which point we work on
after which point doesn't matter
and this updating and calculation scheme is called
Jacoby iteration
the downside of this is it requires lots of iterations
again 024
in this case of a rather course 32 by 32 point crits
luckily there are a few other algorithms which
can speed up its calculation
one of those algorithms being the gossile literation
in Gossidal in contrast to the
to cobiatoration the execution order
so
which point would calculate after which point matters
so let me switch on the point numbers here
and zoom in a bit
what I wanna do for gaocidal is for the interior point
so anything that's not a boundary
start at the lowest point number and calculate its
new value by averaging out the surrounding neighbors
and then taking the next point here
again averaging out the surrounding neighbors
and calculating these values sequentially
the dance to this is that for each point
another one of these neighboring points has
already been updated
so the new values will propagate
faster through our grid
the downside of this of course being its sequential
so I can run it in parallel on all points
let's set this up
and in this case I'll need a few other attributes
so I'm going to need an attribute storing the group
so just if a point is a boundary or not
so let's just wire in a new atrib
create up here
between the group 1 where we create the boundary group
that's this one here
and the group 2
where we create
the group where we set the initial values on
that edge so up here I'll add an atrib create
wire this in here and scrub
only for the boundary group I want to create
that attributor shall call
boundary by default is set to 0
but on this group it is set to 1
and it can be an integer because it's either 01
so after we set this boundary attribute
let's just move all this out of the way
maybe group this call this Jacobi
and further build our gossile
set up after we initialize the points neighbors
and the neighbor numbers
so I want to use an atribrangle here
which I'll set to run over detail
so execute only once
I'm going to wire this up
after we found each points neighbors
and here I want to create a 4 loop which runs
sequentially overall the points
let's call the counting variable PT
for point start at 0 and wild PT
is smaller than the total number of points
that's now PT
let's keep increasing PT by 1 each iterations
okay that's working so far
and in here
the first thing I want to do is create an array
stores a list of
which group is in the boundary group and which isn't
so basically what I wanna do for each of those points
throwing through them sequentially step by step
I wanna read in their boundary attribute
that's the one we just created up here with this
attrib create
and just attach it to a list just add it to a list
for that we'll need a list first
so going to be an integer list let's call it boundaries
and then in our loop here
the first thing I want to do is look at that
boundary attribute
let's call it and boundary
like so it should be equal
to the points boundary
attributed point
coming in through our first input slot
attribute is called boundary
and we're gonna look it up on the point PT again
going through all points sequentially
and then we want to attach
or add this value to the list by pushing it
onto the list boundaries
and we want to push the n bound
re attribute onto that list
and finally let's store this out as a detail attribute
integer list let's call it boundaries
equals to the boundaries list we just created
like so of course I'm missing an at here
and now everything executes fine
okay let's have a look at the Geo spreadsheet
in the detail view here
we can see our first 32 points being in the boundary
and then
following points which are not in the boundary
so I assume this is working so far
let's call this one store boundaries
to list next
I want to do the very same thing to the functions
values to those f values that I created
so basically just copy and paste this node here
and you'll see in a minute why I'm doing this
in two nodes not in one
let's call this stall F2 list
like so and in here I want to change that name
boundaries to F's
Wanna change the output name to F's
change this one to F's and
this one as well
let's call the one we're looking up and F
push it to the list
and instead of the boundary attribute
we're gonna look at the F attribute
which is gonna be a float
float is spells like that
and also the list should be a list of floats like so
okay where do we get wrong here
here I'm trying to store an
integer list when it clearly is a float list
so nothing different than in the previous note
just going through all points sequentially
looking up the F value
pushing it to the list of FS
and then storing that list out as a detail attribute
let's again have a look at the Geo spreadsheed
so down here for the first few points
I'm getting a value of 1
then all zeros
okay that's fine two
back to our note graph
and now it's finally
time to start writing our gossidal update
so again let's use an Attra wrangle
which one set to run only once over detail
and from one of those previous attra brangles here
let's just copy this loop here
paste it down here and delete what's in the loop body
so we just have this loop which runs
over all points sequentially
let's call this note here gossidal update
or Goss on the score Zylo
like so so the first problem we have to solve is
we now have these
values of our points function
of our F values stored in a list
and we haven't stored
any attribute or any list containing
the neighbors of each point
so how do we get those
well let's just switch on the point numbers
have a look up here
maybe to this point here with the number 33
and assuming we stay in that configuration here
where we have a grid which is square
and has the same resolution along its x and Z axis
we can see that it's rather easy to calculate
the point numbers
or list indices of our neighboring points here by
for those two adjacent points here
just subtracting 1
or adding 1 to our current points number
and for these above and below
by subtracting or adding the resolution
remember that is 32
to get those individual point numbers here
to do this first we need to get the grid's resolution
and assuming it's the uniform resolution along the
X and the Z axis
we can just get by with looking up one resolution
so let's call this one RES X
and let's use an integer slider
call it grid resolution
create that slider here
and down here where the slider has been created
I want to link this to the grid's resolution
by going to the grid here
and in this case at the uniform
it doesn't matter if it's rose or columns
it's gonna right click here
and columns go to copy parameter
go down to our gossidal node here
and integrates resolution slider here
I'll right click and paste this as a relative reference
so now we've got 32 in here
also let's load those two lists
the boundaries and the FS list we created here
boundaries and effects
yes so let's first load the boundaries list
let's call it bounds
that should be equal to the list
counteries we stored as a detail attribute
let's also load the values which is a float list
it's called them F's
and they should be equal to the F's list
we just store it
okay back to our senior
let's maybe zoom out again a bit like this
disable point display
and in our loop
the first thing we want to check is if
a point that we're iterating over
is part of the boundary group
and if it is so we want to skip it
so in here we want to first check if
the bounds value at index or list entry PT
is smaller than 1
so if it's 0 and only then
we want to do this
so again we want to look up each points
neighbors values and calculate their average
so let's in this case just go through it manually
yes you could build a loop
but for those four neighboring points
which we have here
I'll stick to manual
so I'm going to specify a flowed variable
which I'll call a val 1
and should be equal to what's stored in our list
for the F's and in here
I want to take my current points number
and from its subtract 1
so that's gonna be my first value
let's just copy that one end paste it down here
and for value 2
I want to just take my point number and to it at 1
let's do it like that
then for the points above and below
remember I want to take
maybe call this about 3
Wanna take my current's point number and to it
add our resolution that we stored
so in our case 32 and let's just copy this line here
paste it down here
and finally for the point below
I wanna subtract this resolution value from it
and store that value in value 4
now let's create the average of that
and that should be a float
and should be equal to val 1 plus val 2 plus val 3
plus val 4
let's put that into parentheses and multiply it by 0.25
which is the same as dividing it by 4
like so and finally
let's set our current points new value
to be this average
so current points value is the same as this average
that's it finally after we've done that
gone through all those points sequentially
let's store out that list again
write it back to a detail attribute
okay not much to see it
maybe
see if we can already see some changed values here
yes you can see there's something going on here
so let's take those values
and write them back from that detail list
onto the individual points using a point wrangle
which goes down here
let's call this one right back F's 2 points
and in here I want to first look up that list of F's
in our case that is a detail attribute
from a current Geo stream which is called F's
and then I want to write back the respective flu
onto my currently active point by just looking up
my currently actives point numbers value
in my list here
like that
and you can see there's something happened again
let's just take this visualization here
where we just offset our grid on why
depending on the value
and we can see after one step we have this
let's save this and
take all
these three notes here and put them in a loop as well
and see how many iterations we need
until we arrive at a decent solution
again using a 4 loop with feedback
gonna cut this
and take the first note that goes
over the store hefts to list
and the last note goes between our visualizer
and the right back
okay 10 iterations
already looking promising
let's maybe go with 32 64
98 to 56 5 12 and 24
and I think after 5 12
there's nothing more happening
or not much more happening
so with that sequential update
we only need to run this loop 512 iterations
versus the
thousand and twenty four of our Jacobi iteration
and let's call this one here
let's put it in a group and call this one gaocidal
alright you might say less iterations
but we can't run
this calcidal update here in parallel
as we could with the jacobiaterations
we doing this sequentially point after point
by running this over
as a detail rangle here
so the speed up might not be that significant
and I would agree with you
however there's one further optimization
which you can do to Gossidal
which will help us
speed up this one
and it's called successive over relaxation or Sor
let's just copy this whole gouseile group and
paste it over
and maybe call this one Sor
what I want to change in here
is only in the gouseile update here
let's call this one gouseile Sor
so instead of blindly calculating
the average value of a point surrounding value
let's instead calculate a correction factor
so the difference that average would make
to my points current value
and for that
just have to look up my points current value
I'm gonna call this Val 0
and it's equal to my current points value
and then after I calculate this average down here
let's calculate the difference
between the Savage value and my points current value
going to take the
average we just calculated and from its subtract val 0
like so and then to just arrive at the same result
what I have to do down here is
not set my current points value as the average
but as its previous value
plus the difference we just calculated like so
and this doesn't make any difference
so
why did we bother doing all this complication in here
well now what we can do is
multiply this difference with a multiplier in order to
over correct and thus taken
into account that not all points in the beginning have
already settled in a sensible position
so multiplying this with a value
between 1 and 2 should result in a quicker correction
let's try a factor of 1.95
which is kind of optimal for this problem
and a good starting value
let's just do that
and you can see our end result didn't change much
however what I can do now
is try decreasing our iterations here
so let's go to 56
you can see still holds up let's decrease it further 2
128 still holds up
to create further to 64 iterations
and now we can see some slight change
so let's try 96
okay a bit of change 128
so 96 iterations looked really decent
now I can see
we massively sped up this calculation by
only requiring 96
iterations of that
sor algorithm
compared to the
224 iterations of our Jacobi and all we do here
is calculate the difference of the
new value that we're calculating towards the old value
and this difference
and each step could be viewed as a correction factor
in this correction factor we just multiply
and thus quote unquote
over correct our solution
to take into account that
not all surrounding points have reached there
final state but still you might say
these futurations are really nice
however for some mesh configurations
this sequentially execution of our update up here
remember this is still a detailer angle
in which we run
overall points sequentially in that loop
so for some configurations large meshes etc etc
this could still in the gate the speed up that we have
by only having this to run
fewer iterations than the Jacobi update
can't we do selective
over relaxation with Jacobi as well
and the answer is by default and no
this is gonna result in errors in your solution
as we're not sequentially updating each
individual point
but there is a way to program around it
and that's what we're gonna cover
in Part 2 of this series
on the persona
and LA Plaza equation
so until then
I hope you had fun with this quite mathy topic
but I hope you'll stay with me
as I think there are some quite
rewarding setups and problems which we can solve
using a LaPlace or poison solver
so see you in Part 2 and as always it is cheers
and goodbye