Free Form Programmable Deformation: Elementary To Finite Element Method (FEM)

Anything you see or sense, anything our machines and sensors scan, are indeed encountering Deformations of many Geometrical structures in nature.?

Simple example is your inner ears; you can hear sounds because the inner ear is highly deformable! The very actions of deformations produce the necessary signals to the brain to hear the outside sounds.?In turn those outside sounds are as well fantastic deformations of the air!

The Rhodopsin molecules in your retina’s sensor cells are deformed upon the impact of a photon i.e. they simply bend! And the consequent quantum mechanical waves serve as signals to the brain reporting seeing something.?

Click on image for more details

Without Deformations we will never be cognizant of the outer world!

Given the Inner ear deformations and the Rhodopsin deformations in Retina we now understand all communications in our universe are direct consequences of Deformation of geometrical structures.

Any animations, any gaming actions are computed by deforming some discretized geometrical objects inside computers.

Click on the image to view the video

1st Order Deformation

The simplest form of algebraic Deformation for parameter 'a' and the increment constant ??a:

a??? ← a? + ??a   n = 0,1,2, ...         

The above reads: if you have the value a? then add the summand ??a and assign the total to a??? .

Example: Assume a = 3 and ??a = 0.236, the following is a sequence defining a's Deformation:

a = 3 and ??a = 0.236

3.000 ← 3 + 0×0.236     0-iteration

3.236 ← 3 + 1×0.236     1-iteration

3.472 ← 3 + 2×0.236     2-iteration

3.708 ← 3 + 3×0.236     3-iteration

...

a??? ← a? + n×??a        n-iteration

or  

a??? ← a? + ??a          n-iteration computes a???        

Example: In Unity this is how 'time' parameter for an animation is incremented:

float time;

void Update( )
{
   time += Time.deltaTime;
   //or
   time = time + Time.deltaTime;
}

// assume a? = 0 is in a??? ← a? + n×??a
// a? ← n×??a

float time=0f; // time begins at 0
int n=-1;      // start with -1 so the first n++ evaluates to 0

void Update( )
{
   n++;
   time = n*Time.deltaTime;
}        

Mass-Spring System and Soft Body Physics

Hooks law's equation for mass-spring is

force = k ??x

where the ??x is the displacement of distance between the two ends of the spring from the rest length.

Q: Where did that formula came from?

A: k ??x is a truncated Taylor Series Expansion of the force function:

Taylor Series Expansion of function f around 0
in interval ??x

f(??x) = c???x? + c???x1 + c???x2 + ... ≈ c???x1 <-- power 1, 1st Order

f(??x=0) = 0 ? c? = 0 which is valid for Spring at rest length ??x=0        

Click on the image to view the video

The grid below represents a matrix/network of orderly connected mass-springs in such a fashion that the initial pulse/force is caused to propagate through the network because of the 1st order Hooks law which induces displacements, eventually, in all mass-springs in the network.

Mass-spring network at rest lengths thus no displacements:

Action of an impulse force causes a local displacement of a few mass-springs:

As time passes the mass-spring system evolves globally meaning all masses displaced in a rather wave-like pattern:

1st Order Deformation by Hook's Law

void Update()
{
   //Loop here to go through all vertices in the mesh
   //??x 
   Vector3 displacement = displacedVertices[i] -
                          originalVertices[i];

   // F=m*??v/??t = ??v/??t  assume m the mass=1
   // Substitute:  k*??x = ??v/??t solve for ??v
   // ??v = k*??x*??t     <--- 1st order Deformation
   // springConstant=k
   // Time.deltaTime=??t
   // displacement = ??x 
   // 0 damper force
   // - or + depending the spring configurations in space
   velocity -= displacement * springConstant * Time.deltaTime;
   ...

}        

2nd Order Deformation and Finite Element Methods

We show how using Free Form Programming Language we are able to solve complicated numerical partial differential equation (PDE) of Standing Waves in 3D dimensions.

u(t, x, y) indicates the z value at time t in (x, y) coordinate.

u(0, x, y) = 0 assuming that at time t=0 all z for all (x, y) points is 0 or all mesh points are lying on a flat x-y plane.

?1u(t=0, x, y) = 0 assuming that at time t=0 all z velocities are 0 at all (x, y) points.

?1u(t) the first derivative with regards to t variable in function u() or z coordinate for the mesh given fixed x and y; computes the displacement of z coordinate u(t+??t)-u(t).

?2u(t) the 2nd derivative with regards to t variable in function u() or z coordinate for the mesh given fixed x and y; computes the change in the displacement of z coordinate by passage of time. In this case we are comparing z coordinates at:

(t, x, y), (t+??t, x y) and (t+2??t, x, y).

You can imagine 2nd derivative with regards to t like unto acceleration and 1st derivative like velocity. At a fixed (x, y) on the mesh computes the acceleration and velocity for z coordinate.

?1u(x) the first derivative with regards to x variable in function u() or z coordinate for the mesh given fixed t and y; computes the displacement of z coordinate u(x+??x)-u(x) where the z at x+??x is compared to z at x, given fixed t and y.

?2u(x) the 2nd derivative with regards to x variable in function u() or z coordinate for the mesh given fixed t and y; computes the change in the displacement of z coordinate where the change in z at x+??x is compared to change in z at x, given fixed t and y. In this case we are comparing z coordinates at:

(x, y), (x+??x, y) and (x+2??x, y).

Same explanation for ?2u(y) given fixed t and x.

The Standing Wave equation is defined as:

?2u(t) - (?2u(x) + ?2u(y)) = 0

It is not difficult to build a model to interpret the 2nd derivatives along x and y axes as four Tangential Tensile vectors; 2nd derivative along the time axis as motion acceleration or again vertical Force F = ma, the wave equation is then rewritten:

??x??y: area of the infinitesimal rectangle

?: area mass density of the membrane

T: Tension force per unit length of the membrane

See Also: 2D Standing Wave Differential Equation, full derivation + 3D mesh geometry of the resulting wave

Since we assumed earlier that the velocity at time 0 being 0 then mesh points only traverse space along the z axis hence the Standing Waves.

Standing waves tend to arise in circumstance where reflection is assured.

Notice the reflections below in a sample Free Form programmable FEM!

Click on the image to view the video

Free Form program which creates and solves and evaluates the Partial Differential Equation (PDE) with Finite Element Method as Method:

region1's type = disk;
region1's radius = 1;
region1's position =(0,0,0);

region2's  type = rectangle
region2's  min point = (-1/2,-1/2);
region2's  max point =(-1/8, -1/8);

region3's type = disk;
region3's radius = 1/7;
region3's position =(1/4,-1/4);

surface waves += finite element with 
equation set to u''(t)-u''(x)-u''(y) = 0 with 

boundary condition set to u=0 with

initial condition = e^(-5(x-0.2)^2+y^2) with

region = region1 - region2 - region3 with 

derivative initial condition set to u'(t)=0 with 

discretization set to method of lines assuming 

max cell measure = 0.1 with 

integration order = 3;        

region1's type: equivalent of region1.type

surface waves: name for the PDE

equation: LHS=RHS form of expression called Equational Expression. u'(x) is shorthand for the first derivative in terms of x, u''(x) second derivative in terms of x. This is an inaccurate syntax. The correct and complete Linear Syntax should be:

derivative(2, 0, 0)(u)(t, x, y)-derivative(0, 2, 0)(u)(t, x, y)-derivative(0, 0, 2)(u)(t, x, y)

or

d(2, 0, 0)(u)(t, x, y)-d(0, 2, 0)(u)(t, x, y)-d(0, 0, 2)(u)(t, x, y)
        

boundary condition: also known as Dirichlet Condition which defines the portions of the boundaries of the region/surface and their shape. u=0 which is actually u(t,x,y)=0 dictates for any point of the boundary of the region to evaluate to z=0 or stationary and flat.

derivative initial condition: the initial values of derivatives of the velocities of the boundary points which is set to u'(t) = or remain with 0 speed and still.

region = region1 - region2 - region3: Subtract "-" operator defines exclusion or cut-out of the LHS from RHS namely cut out or exclude both region2 and region3 from region1.

discretization: Method of Lines is a procedure for discretizing the geometric region for the PDE which is optimized for obtaining fast and accurate solutions.

max cell measure = 0.1: a measure to indicate how tiny the individual cells/polygons of the discretization should be.

integration order: number indicating how many Taylor Series terms used as the defined accuracy requires.

Click on the image to view the video

ASCII JSON Export of FEM Solutions

Per individual mesh vertex, polynomials in t time parameter and an accompanying? boolean expression in t define the procedural animation which was obtained from the solutions of the FEM PDE.

We need to store the output of the FEM solutions and the animations of the membranes into ASCII/JSON formats in order to export into third party systems e.g. Blender or Unity.?

However, there was another reason why we had? to do that: if we stored the FEM meshes obtained from the deformations e.g. at the very minimum mesh deformation we end up with 10s to 100s of megabytes of storage!?

The solution presented itself in Piecewise Polynomial functions which stored only tiny polynomial coefficients such that less than 1 megabyte was consumed to store the symbolic deformations of the FEM meshes solutions.

//A Piece of Piecewise polynomial functions
// if the time is between 0 and 0.0012 seconds this is the function
// for deformation. Thus avoided to store then entire deformed mesh 

{
    0.0160988 -0.000480506 t+0.808115 t^2-815.538 t^3+689721.0 t^4,

    0 <= t <= 0.0012
}        

//Array of symbolic solutions in time per FEM mesh vertex

[
// {polynomial in t, boolean in t} 
{0.0160988 -0.000480506 t+0.808115 t^2-815.538 t^3+689721. t^4,
0<=t<=0.0012},
{0.0160988 -0.000368705 t+0.155944 t^2+581.97 t^3-241951. t^4,
0.0001<=t<=0.002},
{0.0160989 -0.000572152 t+0.775131 t^2-136.729 t^3+34472. t^4,
0.0002<=t<=0.0036},
{0.0160988 -0.000419748 t+0.519041 t^2+35.0385 t^3-4566.07 t^4,
0.0004<=t<=0.0052},
{0.0160989 -0.000508097 t+0.605532 t^2+1.08541 t^3-98.5637 t^4,
0.0008<=t<=0.0068},
{0.0160989 -0.000510342 t+0.606898 t^2+0.752884 t^3-70.8536 t^4,
0.0012<=t<=0.0084},
{0.0160989 -0.000514947 t+0.608675 t^2+0.467957 t^3-54.6645 t^4,
0.002<=t<=0.01},
{0.0160989 -0.000549242 t+0.617806 t^2-0.577638 t^3-11.0981 t^4,
0.0036<=t<=0.0120922}

...

{-117.158+477.582 t-724.441 t^2+485.22 t^3-121.37 t^4,
0.964895<=t<=0.980147},
{-114.166+465.275 t-705.458 t^2+472.207 t^3-118.025 t^4,
0.967946<=t<=0.983198},
{-107.749+438.964 t-665.004 t^2+444.562 t^3-110.941 t^4,
0.970996<=t<=0.986248},
{-102.664+418.182 t-633.148 t^2+422.861 t^3-105.397 t^4,
0.974046<=t<=0.989299},
{-90.5066+368.643 t-557.453 t^2+371.455 t^3-92.3054 t^4,
0.977097<=t<=0.992349},
{-82.8532+337.554 t-510.095 t^2+339.394 t^3-84.1657 t^4,
0.980147<=t<=0.9954},
{-63.3309+258.497 t-390.041 t^2+258.367 t^3-63.6582 t^4,
0.983198<=t<=0.99845},
{-54.4193+222.521 t-335.576 t^2+221.72 t^3-54.4117 t^4,
0.986248<=t<=1.}
...
]        

There will be many many more articles and newsletters on these new Finite Element Methods for Free Form programmable deformations of meshes.

To Be Continued ...

Next Issue: 29th of April 2024