Back to Blog

Realtime Finite Elements Method

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

JerryIJune 20, 2024
differential equationsmodelling3D

Download Notebook

Get the complete code

Initialize mesh

50x50 equally distributed points on 2D plane

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

(*BB[*)(* initial perturbation *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
field2[[25,25]] = 0.05;

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

%3Cspan%20style%3D%22color%3Ared%22%3ERun%20the%20cell%20below%20after%3C%2Fspan%3E evaluation of the initialization cells

Graphics3D[{SpotLight[Blue, {2.4463, 60.1604, 9.3104}, "Intensity"->1000], SpotLight[Red, {48.1378, 47.8679, 12.4545}, "Intensity"->1000], 
  Roughness[0],
  GraphicsComplex[vertices // Offload, Polygon[triangles]],
  EventHandler[AnimationFrameListener[vertices // Offload], 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]];    
  ]]
}, ImageSize->500, Lighting->None, Background->Black]
(*VB[*)(FrontEndRef["1ab83f94-85ab-4a64-aafb-f486cab8adda"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKGyYmWRinWZroWpgmJumaJJqZ6CYmpiXppplYmCUD5RJTUhIBjrwWvg=="*)(*]VB*)

Try to drag it, it is alive