Skip to main content

Realtime Finite Elements Method

· 2 min read
Kirill Vasin

Here we will solve simple 2D wave-equation and visualize it in realtime using polygons

Download original notebook

Initialize mesh​

50x50 equally distributed points on 2D plane

field0 = Table[0., {i,50}, {j,50}];
field1 = field0;
field2 = field0;

(* initial perturbation *)
field2[[25,25]] = 0.01;

lattice = Table[{i,j}, {i, 50}, {j,50}];
mesh = Flatten[lattice, 1];

makeVertices = Compile[{{f, _Real, 2}, {pairs, _Integer, 2}, {scale, _Real}},
  Join[#, {scale f[[#[[1]], #[[2]]]]}] &/@ pairs
];

vertices = NumericArray[makeVertices[field0, mesh, 300]];

Triangulation​

This is necessary for plotting the data using 3D polygons. That we keep indexes of all polygons the same, but will update vertices

Needs["ComputationalGeometry`"];
triangles2[points_] := Module[{tr, triples},
  tr = DelaunayTriangulation[points];
  triples = Flatten[Function[{v, list},
      Switch[Length[list],
        (* account for nodes with connectivity 2 or less *)
        1, {},
        2, {Flatten[{v, list}]}, 
        _, {v, ##} & @@@ Partition[list, 2, 1, {1, 1}]
      ]
    ] @@@ tr, 1];
  Cases[GatherBy[triples, Sort], a_ /; Length[a] == 3 :> a[[1]]]]

triangles = triangles2[mesh];

Let us check the result

ListPlot[mesh, PlotStyle->{Red, PointSize[0.01]}, Prolog->{Line[mesh[[#]]] &/@ triangles}]
(*VB[*)(FrontEndRef["51f56b92-40e0-47fb-8b1a-9df9b700dd1d"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKmxqmmZolWRrpmhikGuiamKcl6VokGSbqWqakWSaZGxikpBimAACAORXR"*)(*]VB*)

Define our equation​

This is discretization of a wave-equation in a form

∂t,t2f−α2∇2f=0\partial^2_{t,t} f - \alpha^2 \nabla^2 f = 0

we follows the following tutorial and rewrite it in a form

update = Compile[{{i, _Integer}, {j, _Integer}, {f, _Real, 2}, {p, _Real, 2}}, 
  0.001 (If[i-1 > 0, f[[i-1, j]], 0] + If[i+1 < 51, f[[i+1, j]], 0] + If[j+1<51, f[[i, j+1]], 0] + If[j-1 > 0, f[[i, j-1]], 0] - 4 f[[i,j]])
+ 2 f[[i,j]] - p[[i,j]]
];

Update function calculated a new value at vertex {i,j}\{i,j\}.

Realtime modelling​

For displaying we use simple Polygon and GraphicsComplex

Run the cell below after
evaluation of the initialization cells

Graphics3D[{Shadows[True], SpotLight[Blue, {2.4463, 60.1604, 9.3104}, 1.5, 3000], SpotLight[Red, {48.1378, 47.8679, 12.4545}, 1.5, 3000], 
  Roughness[0],
  GraphicsComplex[vertices // Offload, Polygon[triangles]],
  AnimationFrameListener[vertices // Offload, "Event" -> "frame calc"]
}, ImageSize->500]
(*VB[*)(FrontEndRef["93deb6c9-ae46-42de-b1aa-0dd3fdbb85eb"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWxqnpCaZJVvqJqaamOmaGKWk6iYZJibqGqSkGKelJCVZmKYmAQCW7Bbt"*)(*]VB*)

try to pan and zoom it ! 😉

Now we wrap our calculations into a animation cycle

Evaluate this if animation did not start

(*BB[*)(*reset again*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
field0 = Table[0., {i,50}, {j,50}];
field1 = field0;
field2 = field0;

(* initial perturbation *)
field2[[25,25]] = 0.01;

EventHandler["frame calc", Function[Null, 

  With[{i = #[[1]], j = #[[2]]}, 
    field0[[i,j]] = update[i,j, field1, field2];
  ]&/@ RandomSample[mesh];
  
  field2 = field1;
  field1 = field0;
  
  vertices = NumericArray[makeVertices[field0, mesh, 300]];

]];

(*BB[*)(*kickstarter*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
EventFire["frame calc", True];

To stop it

EventRemove["frame calc"]